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Wide-angle  elastic  wave  one-way  propagation  in  heterogeneous 
media  and  an  elastic  wave  complex-screen  method 

Ru-Shan  Wu 

Institute  of  Tectonics,  University  of  California,  Santa  Cruz 


In  this  paper  a  system  of  equations  for  wide-angle  one-way  elastic  wave  propagation  in  arbitrarily  hetero¬ 
geneous  media  is  formulated  in  both  the  space  and  wavenumber  domains  using  elastic  Rayleigh  integrals  and 
local  elastic  Bom  scattering  theory.  The  wavenumber  domain  formulation  leads  to  compact  solutions  to 
one-way  propagation  and  scattering  problems.  It  is  shown  that  wide-angle  scattering  in  heterogeneous  elastic 
media  cannot  be  formulated  as  passage  through  regular  phase-screens,  since  the  interaction  between  the  in¬ 
cident  wavefield  and  the  .heterogeneities  is  not  local  in  both  the  space  domain  and  the  wavenumber  domain. 
Our  more  generally  valid  formulation  is  called  the  "thin-slab”  formulation.  After  applying  the  small-angle 
approximation,  the  thirt-slab  effect  degenerates  to  that  of  an  elastic  complex-screen  (or  "generalized  phase- 
screen").  Compared  with  scalar  phase-screen,  the  elastic  complex-screen  has  the  following  features.  (1)  For 
P-P  scattering  and  S-S  in-plane  scattering,  the  elastic  complex-screen  acts  as  two  separate  scalar  phase- 
screens  for  P  and  S  waves  respectively.  The  phase  distortions  are  determined  by  the  P  and  S  wave  velocity 
perturbations  respectively.  (2)  For  P-S  and  S-P  conversions,  the  screen  is  no  longer  a  pure  phase-screen  and 
becomes  complex  (with  both  phase  and  amplitude  terms);  both  conversions  are  determined  by  the  shear  wave 
velocity  perturbation  and  the  shear  modulus  perturbation.  For  Poisson  solids  the  S  wave  velocity  perturba¬ 
tion  plays  a  major  role.  In  the  special  case  of  OCq  =  2p0,  S  wave  velocity  perturbation  becomes  the  only  fac¬ 
tor  for  both  conversions.  (3)  For  the  cross-coupling  between  in-plane  S  waves  and  off-plane  S  waves,  only 
the  shear  modulus  perturbation  8jl  has  influence  in  the  thin-slab  formulation.  For  thexomplex-screen  method 
the  cross-coupling  term  is  neglected  because  it  is  a  higher  order  small  quantity  for  small-angle  scattering. 
Relative  to  prior  derivations  of  vector  phase-screen  method,  our  method  can  correctly  treat  the  conversion 
between  P  and  S  waves  and  the  cross-coupling  between  differently  polarized  S  waves.  A  comparison  with 
solutions  from  three-dimensional  finite  difference  and  exact  solutions  using  eigenfunction  expansion  is  made 
for  two  special  cases.  One  is  for  a  solid  sphere  with  only  P  velocity  perturbation;  the  other  is  with  only  S 
velocity  perturbation.  The  Elastic  complex-screen  method  generally  agrees  well  with  the  three-dimensional 
finite  difference  method  and  the  exact  solutions.  In  the  limiting  case  of  scalar  waves,  the  derivation  in  this 
paper  leads  to  a  more  generally  valid  new  method,  namely,  a  scalar  thin-slab  method.  When  making  the 
small-angle  approximation  to  the  interaction  term  while  keeping  the  propagation  term  unchanged,  the  thin- 
slab  method  approaches  the  currently  available  scalar  wide-angle  phase-screen  method. 


1.  Introduction 

The  one-way  wave  equation  has  been  widely  used  for  calcu¬ 
lating  acoustic  wavefields,  especially  in  ocean  acoustics.  The 
one-way  wave  equation  is  a  marching  algorithm  similar  to  a 
Markov  process,  in  which  the  calculation  of  the  next  wavefront 
requires  only  knowledge  of  the  current  wavefront.  Compared 
with  the  two-way  full  wave  equation,  in  which  at  any  moment 
the  wavefield  is  related  to  the  wavefield  of  the  previous  time 
step  for  the  whole  space,  the  one-way  equation  approach  has 
the  advantages  of  fast  computation  speed  and  less  storage  re¬ 
quirement,  and  at  present  is  the  only  method  to  calculate  the 
long-range  acoustic  wave  propagation  in  the  range-dependent 
ocean  environment  [e.g.,  Tappert ,  1977;  Platte' and  Tappert , 
1975;  Collins ,  1989,  1990].  In  exploration  geophysics,  it  has 
been  used  for  both  migration  and  modeling  [e.g.,  Stoffa  et  aL, 
1990;  Wu  and  Huang ,  1992]. 

The  generalization  of  the  one-way  wave  equation  to  elastic 
waves  has  great  potential  for  applications  in  many  disciplines, 
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such  as  modeling  and  inversion  in  seismic  exploration,  and 
long  range  seismo-acoustic  wave  propagation  in  the  ocean  with 
the  sedimentary-bottom  and  ice-cap  interactions.  Modeling 
long-range  crustal  wave  propagation  in  heterogeneous  wave 
guides  to  study  the  propagation  and  scattering  of  regional 
phases  such  as  Lg,  which  is  an  extensively  used  wave  group 
for  yield  estimation  of  underground  nuclear  explosions,  presents 
another  need  for  such  a  technique.  Especially  for  three- 
dimensional  (3D)  cases,  the  full  two-way  elastic  wave  finite 
difference  method  (for  a  review,  see  Frankel  [1989]  is  not  real¬ 
istic  for  the  study  of  long-range  propagation  and  scattering 
effects. 

The  traditional  way  of  deriving  the  one-way  parabolic  equa¬ 
tion  for  scalar  waves  is  to  assume  a  solution  in  the  form  of 
p  -  A  exp (ikxx),  where  kx  is  the  principal  wavenumber  in  the 
propagation  direction  and  A  is  the  slowly  varying  amplitude 
function.  One  substitutes  this  trial  solution  into  the  two-way 
full  wave  equation.  Taking  the  small-angle  approximation  by 
dropping  higher  order  derivatives  in  the  propagation  direction 
results  in  a  parabolic  equation  for  the  amplitude  function  (see 
Tappert ,  [1977].  This  approach  became  very  popular  in  acous¬ 
tics.  The  generalization  from  scalar  wave  to  elastic  wave  para¬ 
bolic  equations  was  first  introduced  by  Landers  and  Claerboul 
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[1972]  for  a  two-dimensional  (2D)  case  along  similar  lines.  In 
their  derivation,  the  dilatation  and  the  single  nonzero  com¬ 
ponent  of  the  rotation  of  the  elastic  wavefield  were  used  as 
variables,  and  a  higher  order  approximation  was  made  to  ex¬ 
tend  the  valid  range  to  wider  propagation  angles.  Hudson 
[1980]  used  the  displacement  fields  as  variables  in  his  deriva¬ 
tion  but  followed  the  same  single  principal  wavenumber  ap¬ 
proach.  His  derivation  is  for  the  3D  case  but  the  approximation 
was  made  only  for  small-angle  propagation  (parabolic  approxi¬ 
mation).  Now  we  know  that  this  single  principal  wavenumber 
approach  is  only  good  for  the  common-mode  scattering  prob¬ 
lem  and  is  intrinsically  not  appropriate  to  deal  with  wave 
conversion  problems.  In  the  case  of  elastic  media,  there  exist 
two  wave  speeds,  the  compressional  and  the  shear  wave  speeds. 
Therefore,  two  separate  parabolic  equations  must  be  obtained 
for  both  the  P  and  S  waves.  Although  each  equation  has  a 
coupling  term,  this  pseudo*-con  version  term  may  not  correctly 
represent  the  real  physical  conversion,  because  the  neglected 
terms  during  the  process  of  making  the  parabolic  approximation 
are  not  necessarily  small  compared  with  the  remaining  terms 
due  to  the  large  difference  between  the  P  and  S  wavenumbers. 
It  has  been  shown  that  in  some  cases,  this  approach  can  lead  to 
unreasonably  large  errors  [Wales  and  McCoy ,  1983]. 

Another  approach  is  to  derive  the  one-way  equation  for  lay¬ 
ered  media  (one-dimensional  media).  In  this  case  the  downgo¬ 
ing  and  upgoing  waves  are  decoupled;  therefore  a  one-way  ma¬ 
trix  equation  can  be  obtained  (see:  Ursin ,  1983;  Kennett,  1983; 
Wapenaar  and  Berkhout,  1989;  Frazer ,  1990].  However,  this 
approach  does  not  have  the  capacity  of  dealing  with  3D  propa¬ 
gation  and  scattering  problems.  Some  types  of  one-way  finite 
difference  equations  can  be  also  derived  for  the  case  of  a  lay¬ 
ered  medium,  which  can  be  applied  for  weakly  range-dependent 
propagation  problems  [ Greene ,  1984,  1985;  Wetton  and  Brooke , 
1990;  Collins ,  1989,  1991].  In  this  approach,  the  vertical  dis¬ 
placement  and  the  dilatation  are  used  as  variables.  After  nor¬ 
malization  of  variables  for  cylindrical  spreading,  the  governing 
differential  equations  become  2D  coupled  equations,  which  in 
turn,  are  factored  into  one-way  equations  and  the  square-root 
operator  is  approximated  by  a  higher  order  Pade "approximation, 
resulting  in  a  system  of  stable  elastic  parabolic  equations.  The 
equations  are  solved  by  a  finite  difference  algorithm  and  can 
be  used  for  wide-angle  and  weak  range-dependent  propagation 
problems.  However,  this  formulation  is  not  capable  of  dealing 
with  many  scattering  problems,  such  as  the  SV-SH  cross¬ 
coupling,  strong  range  variations  and  real  3D  scattering  prob¬ 
lems. 

McCoy  [1977]  and  Wales  [1986]  adopted  a  perturbation  ap¬ 
proach  and  derived  a  system  of  elastic  parabolic  equations  for 
the  P  and  S  wave  potentials.  The  two  wave  speeds  are  put  into 
the  formulation  in  the  beginning;  therefore  the  wave  conversion 
problem  can  be  treated  properly.  The  final  result  is  in  the  form 
of  a  marching  finite  difference  algorithm  in  the  space  domain. 
The  method  is  suitable  for  small-angle  forward  propagation 
problems.  Since  the  free  boundary  conditions  for  the  potentials 
are  rather  difficult  to  formulate,  there  are  some  stability  prob¬ 
lems  for  certain  wave  types  [Wales,  1986]. 

Fisk  and  McCartor  [1991],  hereafter  identified  as  FM9I, 
have  introduced  a  vector  phase-screen  method  for  elastic 
waves,  and  made  some  comparisons  with  the  2D  finite 
difference  method  [Fisk  et  ai ,  1992].  They  assumed  that  the 
medium  is  smooth  and  the  gradients  of  parameter  variation  can 
be  neglected,  resulting  in  a  pair  of  decoupled  P  and  S  wave 
equations.  The  scattering  effects  are  simulated  by  passing  the  P 


and  S  waves  through  two  separate  phase-screens  involving  P 
and  S  wave  speed  perturbations,  respectively.  The  P-S 
conversion  is  obtained  by  projecting  the  distorted  P  wavefront 
to  the  corresponding  S  wave  components.  The  same  is  done  for 
the  S-P  conversion.  The  procedure  is  relatively  simple.  How¬ 
ever,  there  exists  some  intrinsic  difficulty  and  internal  incon¬ 
sistency  in  this  approach.  First,  the  basic  assumption  of  decou¬ 
pled  P  and  S  wave  equations  implies  that  the  theory  is  not  ca¬ 
pable  of  dealing  with  the  conversion  problem.  The  theory  can 
only  handle  the  first-order  effects  of  the  wave  front  distortion 
due  to  velocity  perturbations,  such  as  focusing  and  defocusing, 
and  coda  generation  by  forward  scattering.  Formally  projecting 
the  distorted  P  wave  front  into  S  wave  components  may  pro¬ 
duce  nonphysical  conversions.  We  will  discuss  this  in  section 

4.  Second,  there  is  ambiguity  about  the  scattering  effects  of 
perturbations  of  different  parameters.  In  their  method  (FM91), 
P-P  and  P~S  scattering  are  determined  by  P  wave  velocity 
perturbations  5a,  while  SS  and  S-P  scattering  are  deter¬ 
mined  by  S  wave  perturbations  5(3,  leading  to  ambiguity  in 
selecting  the  combination  of  three  perturbations  8p,  5X,  and  8)i 
for  an  isotropic  elastic  medium.  The  P  wave  velocity  perturba¬ 
tion  5a  could  be  an  arbitrary  combination  of  5p,  5X,  Sji.  While 
these  perturbations  can  have  very  different  scattering  effects, 
the  method  they  adopted  could  not  distinguish  between  these 
effects,  which  implies  again  that  the  method  cannot  handle  the 
wave  conversion  properly. 

In  this  paper  I  will  derive  one-way  elastic  wave  equations  for 
P  and  S  displacement  fields  using  the  elastic  wave  Rayleigh  in¬ 
tegrals  and  elastic  Bom  scattering  theory  following  a  perturba¬ 
tion  approach.  First  the  wide-angle  "thin-slab"  formulas  are 
derived  in  space  domain  (section  2)  and  wavenumber  domain 
(section  3).  The  method  can  be  implemented  in  the 
wavenumber  domain  similar  to  a  generalized  split-step  algo¬ 
rithm  for  elastic  waves,  and  applied  to  wide-angle  fonvard 
scattering  problems.  Next,  an  elastic  "complex-screen"  algo¬ 
rithm  is  obtained  by  using  the  small-angle  approximation  (para¬ 
bolic  approximation)  (section  4),  which  can  be  implemented  by 
a  dual-domain  technique  using  fast  Fourier  transform  (FFT), 
which  has  much  faster  speed  (2  to  3  orders  of  magnitude)  than 
the  thin-slab  method.  Computation  speed  is  discussed,  and 
comparison  with  finite  difference  and  exact  solutions  for  an 
elastic  sphere  is  made  through  numerical  examples  (section  4). 
Finally,  the  limiting  case  of  scalar  waves  is  obtained  in  section 

5.  This  approach  leads  to  a  scalar  thin-slab  algorithm,  which  is 
valid  under  more  general  conditions  than  the  currently  available 
"wide-angle"  phase-screen  methods.  When  making  the  small- 
angle  approximation  to  the  interaction  term  while  keeping  the 
propagation  term  unchanged,  the  method  approaches  the  regular 
wide-angle  scalar  phase-screen  [Feit  and  Fleck ,  1978;  Thomson 
and  Chapman ,  1983], 

2.  Space  Domain  Formulation 

General  Media 

The  equation  of  motion  in  a  linear,  heterogeneous  elastic 

p(x)il(x)  =  V-  a(x)  (1) 

where  u  is  the  displacement  vector,  o(x)  is  the  stress  tensor 
(dyadic)  and  p  is  the  density  of  the  medium.  Here  u  is  defined 
as  d2u fdt2.  We  assume  no  body  force  exists  in  the  medium.  We 
know  the  stress-displacement  relation 
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or(x)  =  c(x) :  £(x)  =  ~c :  (Vu  +  u V)  (2) 

where  c  is  the  elastic  constant  tensor  of  the  medium,  £  is  the 
strain  field,  uV  stands  for  the  transpose  of  Vu,  and  stands 
for  double  scalar  product  of  tensors  defined  through 
(ab):(cd)  =  (bc)(a*d).  Equation  (1)  can  then  be  written  as  a 
wave  equation  of  the  displacement  field: 

p(x)ii(x)  =  V*  [yc :  (Vu  +  uV)]  .  (3) 

If  the  parameters  of  the  elastic  medium  and  the  total  wave 
field  can  be  decomposed  as 

P(x)  =  p0  +  Sp(x) 

c(x)  =  Co  +  5c(x) 

u(x)  =?  u°(x)  +  U(x)  (4) 

where  p0  and  Cq  are  the  parameters  of  the  background  medium, 
5p  and  5c  are  the  corresponding  perturbations,  u°  is  the 
incident  field  and  U  is  the  scattered  field,  then  (3)  can  be 
rewritten  as 


p0U-V*  [yCo :  (VU  +  UV)]  =  Q  (5) 

Q  =  -|spu  -  V-  [y 5c :  (Vu  +  uV)]j  =  -  5pii  +  V-  [5c :  £]  (6) 

where  Q  is  the  equivalent  body  force  due  to  scattering.  Note 
that  u  also  satisfies  equation  (5). 

If  the  wave  field  is  known  on  a  closed  surface  5,  then  from 
the  representation  theory  [Aki  and  Richards ,  1980]  we  can 
express  the  solution  of  (5)  as  the  sum  of  a  surface  integral  and 
a  volume  integral 

u(x0  =  <J(x)]*  G(x,x!)-u(x>  [j h •  £(x,X!)]| 

+  ^dV  Q(x)«  G(Xl,x)  (7) 

where  a(x)  is  the  stress  field  on  the  surface,  T(x)  =  n-  a(x)  is 
the  corresponding  traction  field,  u(x)  is  the  displacement  field 
on  the  surface,  n  is  the  outward  normal,  G(x,xj),  ZfoxO,  and 

2-n(x»Xi)  =  n •  X(x,xj)  are  the  Green’s  displacement,  stress,  and 
traction  tensors,  respectively.  This  representation  integral  in 
the  present  form  was  first  introduced  and  discussed  by  Love  in 
1904  (in  a  more  general  form  first  by  Betti  in  1872  [e.g.,  Aki 
and  Richards  1980]),  and  is  called  the  Love  integral  or  the 
elastic- wave  Kirchhoff  integral  [ Pao  and  Varatharajulu,  1976]. 
Note  that  the  Love  integral  (the  surface  integral)  is  from  the 
contribution  of  outside  sources,  i.e.  the  sources  outside  the  sur¬ 
face.  The  known  field  on  the  surface  can  be  generated  by  actual 
sources  outside  the  surface,  or  by  scattering  of  the  hetero¬ 
geneities  outside  the  surface.  The  scattered  field  generated  by 
the  heterogeneities  inside  the  surface  will  have  zero  value  of 
the  surface  integration  and  therefore  does  not  contribute  to  the 
Love  integral.  We  can  then  call  the  part  of  field  represented  by 
the  Love  integral  as  the  primary  field.  Therefore,  (7)  becomes 


u(x,)  =  lAxO  +  UtxO  (g) 


In  equation  (10)  integration  by  parts  has  been  used  so  that  the 
volume  integration  contains  no  gradients  of  the  elastic  constants 
[see  Wu  and  Aki ,  1985;  Wu ,  19896],  Note  that  in  (7)  and  (10) 
the  integration  volume  V  is  over  the  whole  volume  which  con¬ 
tains  all  the  scatterers.  When  performing  integration  by  parts, 
the  surface  is  chosen  to  be  outside  the  heterogeneous  region,  so 
that  the  surface  terms  vanish. 

Now  we  formulate  the  representation  integral  for  the  for¬ 
ward  propagation  problem.  Suppose  the  incident  wave  is  pro¬ 
pagating  along  the  x-axis.  Now  consider  that  the  integration 
surface  in  the  Love  integral  is  composed  of  a  plane  surface  S 
at  x,  which  is  perpendicular  to  the  propagation  direction,  and  a 
hemisphere  with  a  large  radius.  As  seen  in  Figure  1,  we  can 
assume  that  the  major  contribution  to  the  surface  integral  to 
calculate  the  primary  field  between  planes  S  and  is  from  the 
large-aperture  plane  surface  5.  The  contribution  from  the  other 
parts  of  the  closed  surface  can  be  neglected  due  to  the  weak¬ 
ness  and  time  delay  of  the  signal.  Therefore  the  Love  integral 
gives  the  elastic  wave  field  generated  by  the  sources  and  by 
wave  scattering  of  all  the  heterogeneities  on  the  left-hand  side 
of  plane  S. 

Now  let  us  examine  the  contribution  from  volume  hetero¬ 
geneities.  Now  the  volume  integral  accounts  for  the  contribu¬ 
tion  from  all  the  scatterers  on  the  right-hand  side  of  S .  How¬ 
ever,  under  the  one-way  propagation  assumption,  the  backscat- 
tered  waves  are  neglected,  and  therefore  the  contribution  from 
the  scatterers  on  the  right-hand  side  of  S\  on  which  the  obser¬ 
vation  point  Xi  is  situated,  is  excluded.  Furthermore,  if  the  dis¬ 
tance  between  S  and  Si  is  small,  then  the  field  u(x'),  where  x' 
is  situated  between  x  and  Xj,  can  be  approximated  by  u°(x'), 
since  the  modification  of  the  wave  field  due  to  scattering 
between  S  and  S  i  is  small  compared  with  u°(xO.  Therefore,  the 
scattered  field  can  be  approximated  as 


Fig.  1.  Geometry  of  the  derivation.  The  medium  is  sliced  into  thin- 
slabs.  Now  consider  the  scattering  and  propagation  effects  of  the  thin- 
slab  between  plane  S  and  S\.  Suppose  the  wave  field  is  known  on  S. 
The  wave  field  on  S  j  will  be  formulated  by  the  forward  propagation 
approximation. 
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where  V'  means  the  operation  is  taken  with  respect  to  x\  a 
point  inside  the  slab.  In  fact,  this  is  an  integral  of  elastic  Bom 
scattering  [Gubernatis  et  al ,  1977;  Wu  and  Aki ,  1985;  Wu , 
19896].  Here  we  can  consider  (11)  as  a  local  Bom  approxima¬ 
tion. 

In  order  to  calculate  u(xt),  we  need  to  determine  u°(xi)» 
u°(x')  and  £  (x')  from  the  known  field  on  S  using  the  Love 
integral.  It  has  been  shown  that  in  the  case  of  one-way  propa¬ 
gation  from  a  plane  surface,  the  Love  integral  can  be  replaced 
by  an  elastic-wave  Rayleigh  integral  (ERI)  which  is  simpler 
and  computationally  much  more  efficient  [Wu,  1989a; 
Wapenaar  and  Berkhout ,  1989] 

u°(xi)=  -2|^Su(x>Sn(x,x1).  (12) 

This  is  ERI  of  type  II.  Similarly,  we  can  calculate  the  strain 
field 

£°(x')  =  -  2j dS  u(x>  £,'[£„  (x,xO]  ( 1 3) 

where  8X[*  ]  is  the  strain  operator: 

£x[a]  =  IfV.a  +  aVJ  .  (14) 

where  Vx  stands  for  taking  a  gradient  with  respect  to  x . 

Equivalently,  if  the  traction  field  T(x)  is  known  on  S ,  the 
strain  field  can  be  calculated  using  the  ERI  of  type  I: 

£°(x')  =  2 Jc/S  T(x>  E'(x,x')  (15) 

where  E'(x,x')  is  Green’s  strain  tensor,  defined  as 

E'(x,x')  =  £x'[G(x,x')]  =  y  [V'G(x.x')  +  ( V'G(x,x'))2 13]  ( 1 6) 

where  V'G  stands  for  taking  the  gradient  of  G  with  respect  to 
x'  and  (V'G)213,  the  transpose  of  V'G  with  respect  to  the  first 
two  indexes.  The  traction  field  is  related  to  the  strain  field  by 


Assuming  a  homogeneous  background.  Green’s  tensors  in  the 
above  formulas  can  be  written  out  explicitly.  For  details  see 
Appendix  A. 

3.  Wave  Number  Domain  Formulation 

Referring  to  the  geometry  in  Figure  1,  Let  us  calculate  the 
scattered  field  by  the  heterogeneities  within  a  vertical  thin  slab 
of  width  Ax  =  xx-x.  Assume  we  know  the  displacement  field 
on  the  surface  S,  u°(x,xr),  where  xT  -  (y,z)  is  the  2D  position 
vector  in  the  y-z  plane.  For  the  sake  of  simplicity,  we  put 
x  =  0  in  the  following  derivation. 

We  decompose  u°  on  the  surface  S  into  plane  waves 

u°(0,xr)  =  Ua(0,xr)  +  Up(0,xr) 

u°(0,xr)  =  -L-pKr  [u°(Kr)  +  Up°(Kr)kiV  ‘r  .  (22) 
47T  j 

Then  the  displacement  field  u°  everywhere  in  the  slab  can  be 
expressed  as 

u°(x')  =  u«(x')  +  u£(x') 

=  -—UKt  [u°(Kr)  A' x'  +  up°(Kr)eikl1' X']  (23) 

47T  j 

where 

ka  =  -Kt*x  +  =  ya4  +  Kr 

kp  =  -K}ex  +  Kr  =  +  Kr  (24) 

and  the  strain  field 

eV)  =  £„(x')  +  Ep(xO  =  L(V'u0  +  u0V') 

=  -L-LKr  [ika(£au°  +  u°£a)e‘ka  x 
871^  J 

+  z/:p(^pu^  +  u^p)^'kp  *].  (25) 


T(x)  =  n  •  [c(x) :  £(x)] .  (17) 


Isotropic  Media 

In  the  case  of  isotropic  media,  equations  (8),  (12),  (11),  and 
(15)  can  be  further  simplified  to  (in  the  frequency  domain) 

u(x,)  =  U°(xi)  +  U(x0  (18) 

u°(x,)=  -2jdSu(x).ZfI(x,x1)  (19) 

U(x,)  =  co 2(dV  8p(x')u°(x0*  G(x!,x') 


+  \dV  [5>.(x')! £°(x')| I  +  25p(x')£°(x') 


:  V'G(xi.x')  (20) 


with  u°(x')  and  £  (x')  calculated  by  ERI  of  type  II  and  type  I, 
respectively.  For  the  calculation  of  £  (x'),  (15)  can  be  used 
and  the  traction  field  can  be  determined  in  this  case  by 


T(x)  =  n  * 


X(xi  £(x)(  I  +  2p(x)£(x)J 


(21) 


where  I  is  the  unit  tensor  (idemfactor)  and 
|£(x)|  =  I: £(x)  =  V*  u(x). 


Also  we  decompose  the  scattered  field  U(x)  expressed  by  (20) 
on  the  surface  Si  into  plane  waves.  Substituting  (23),  (25),  and 
the  wavenumber  domain  expressions  of  Green’s  tensors  (see  the 
derivation  in  Appendix  B)  into  (20),  we  can  derive  the  expres¬ 
sions  of  the  elastic  wave  Bom  scattering  in  wavenumber 
domain.  As  a  demonstration  of  the  procedure,  we  derive  here 
the  scattered  P  wave  contributed  by  5p  for  an  incident  P 
wave: 


U£p(x)  =  (i>2^dV  8p(x')u°(x>  Ga(x,x') 

=  ^KVpKr|^5p(x'),,k“'X' 


ikl 


'  87t2p0C0Ya 


(u«(Kr)-  k'a)k'ae 


i  k'  *  (x-xO 


47t2  J 


[dK're'^-’rll 


l'y'nAx 


2  Ya 


x-LfdKT  k25p(k'tt:-~-u°(KT)(ka-  k'a)k'a  (26) 
4n  J  Po 

where  Ua(KT)  =  lua(Kr)l  and  5p(ka)  is  the  3D  Fourier 
transform  of  5p(x): 
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5p(ka)  =  | dV  5p(x,xr)e  ‘  7  *r, 


,  kt*  xre  *V 


=  5pQt,Kr)e  v dx  . 


Therefore, 


W  ».2  1 


KV)  = 

..  5p(Ya-Ya.K'r  -Kr)  0/n  „ 

X  Me  (OjIvt" )(Xa*  a  oJX  a  .  (28) 

Po 

/  y  r 

We  recognize  that  the  phase  term  e  07  represents  the  plane 
wave  propagation  between  the  two  surfaces.  Then  the  5p/p 
contribution  to  the  scattered  P  wave  of  the  wavenumber  K'r 
for  an  incident  plane  wave  with  wavenumber  Kr  can  be 
represented  as 

U"(K'r,Kr)  =  -L^--P(k'°'ka)^(Kr)(tV  *'«)£'«  (29) 
2y  a  po 

This  is  the  wavenumber-domain  formula  of  elastic  wave 
scattering  by  the  Bom  approximation  (here  the  local  Bom 
approximation).  It  is  similar  to  that  of  the  far-field  Bom  scatter¬ 
ing  [see  Wu  and  Afci ,  1985;  Wu ,  19896]  except  for  a  weighting 
factor  i/2/a.  However,  in  this  wavenumber-domain  formula¬ 
tion,  no  far-field  approximation  is  made.  The  only  approxima¬ 
tion  is  the  local  Bom  approximation. 

Following  the  same  procedure  for  other  parts  of  the  scattered 
field,  we  have 

U(K'r,Kr)  =  Uw(KV,Kr)  +  Uw(KV,Kr) 

+  Uw(KV,K70  +  Uw(KV,K7-)  (30) 

U"(KV,Kr)  =  -^7-klulk'l^  ^Ml 

2y  a  I  Po 


8^(£)  (£  .  £,  \2  25p(k) 

Jo  +  2|io’1  a  j  Xo  +  2po 


Uw(K'r,K r)  -  J-klu^[ka-kf^  ^ 

2Yp  p  Po 


where  =  !u£(Kr)l,  uj?  =  u$(Kr),  a  =  u£/iuj?l  and 
m  =  M/iMl,  M  =  u£-j£'p(u£’  if'p) 

f  -  m  X^'p  . 

In  (31),  k  =  kJC  -km  is  the  exchange  wavenumber  with  k"\  k*c 
as  the  incident  and  scattering  wavenumbers,  respectively,  and 

K  =  V*a  -K?ex  +  Kr 

kp  =  'Jkf-Kjex+  Kt 

k'a  =  V^'-ATV^+KV 

k'p  =  V^p2  +K'r  •  (32) 

Note  that  m  is  perpendicular  to  the  outgoing  direction  £'p  and 
lies  in  the  plane  determined  by  Up  and  £'p.  Therefore,  the  first 
term  of  the  last  equation  in  (31)  stands  for  the  in-plane  S  wave, 
i.e.  the  non-coupled  S  wave;  while  /  is  the  off-plane  unit  vec¬ 
tor  which  is  perpendicular  to  both  the  outgoing  direction  and 
m.  We  see  the  cross-coupling  (the  depolarization  of  S  waves) 
can  only  be  generated  by  shear  modulus  perturbations  5)1. 

To  obtain  the  scattered  field  at  the  surface  Sx  with  transverse 
wavenumber  K'r>  we  need  to  sum  up  the  contributions  from  all 
the  incident  plane  waves: 


1 

\d  Kre^ 

4k2  • 

)  1 

1 

[ dKTe *** 

4k2' 

)  1 

where  Ax  is  the  distance  between  the  two  surfaces.  To  account 
for  the  scattering  effect  of  the  heterogeneities  between  x  and  x\ 
the  integration  limits  in  (27)  need  to  be  changed  to  x  and  x\. 

The  free  propagation  of  the  incident  field  in  the  wavenumber 
domain  can  be  obtained  from  (19)  by  substituting  (23)  and  the 
plane  wave  decomposition  of  (x,xO  (Appendix  B): 

u°(xi)  =  2\dS  u(x)-  Xx(x,x,) 


-  2(— )(£„•  ^Mi 

ao  Po 


U"  (KV.Kr)  =  -±-kl{ up°.  f'jjf'jML 
2Ta  Po 


•2(— )(£„•  jf'jMl 
«0  Po 


=  2|^|u5(x)*  Zax(x,Xj)  +  up(x)*  Zpx(x,xj) 

=  --t-tiS  [dKT  \dK'T  }  — * 

2n2i  J  rJ  r87t2po)2l 


~e  kV  ^ ' "e ' k°‘  *  X(u°-  x )k'a  +  2p(u°-  k'a)(x •  if'Jif'J 

la 


+  Vie/kV  (v*w^ 

Yp 


Uss (K't,Kt)  =  Jf'p)]M- 

-  .(V  ^'p)(up°  -  £'p(up-  £'„)) + (up0-  £'p)(4  -  fyjfp.  if '„))]  Ml 


X  [(*  •  £ ‘ 'p)up°  +  (up0-  £ 'p)x  -  2(x  •  k‘ 'p)(up°-  'p  | .  (34) 

Knowing  that 


-  l( up0-  *'p)(/-  jfcp)Mi 


pei(k“'k“)'X  =  ei<Y“-y“,X|^ei(Kr' 


=  47i25(Kr  -  K'r)e'<ra  ^a)x 


jdS  /  W  x  =*4jt25(Kr  -  ^  (35) 
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and  noticing  that  (up*  k$)  -  0,  the  integration  over  K'r  can  be 
completed,  resulting  in 


u°(x.) 


1 _ 

87l2pC02 


X|a  +  2p)(C^)Ua°(Kr) 


+  —  e'kp  *V(x*  ^p)up(Kj-) 

Yp 


(36) 


By  inspection,  we  obtain 

uVi,Kr)  =  K(x,KT)eiyaAx  +  upV.K^e'^  .  (37) 


As  expected,  this  is  simply  the  free  propagation  between  the 
two  planes. 

Finally,  the  total  field  at  is  obtained  by  summing  up  the 
primary  field  (37)  and  the  scattered  field  (33): 


it  (xlfKV)  =  e  'a  W(xXt) 


+  -CpKr  [u"(K'r,Kr)  +  Uw(K'r,Kr) 

471  J 


us(x„K'r)  =  e  P  < u$(x ,K't) 


X~t\dKT  [uM(KV,Kr)  +  UM(KV,Kr) 

TTZ  ^ 


(38) 


The  above  procedure  can  be  implemented  iteratively  in  the 
wavenumber  domain,  and  the  propagation  will  march  step  by 
step  in  the  forward  direction.  For  each  new  step,  the  total  field 
calculated  using  (38)  from  the  previous  step  is  taken  as  the 
incident  field,  and  then  the  scattered  and  primary  fields  on  the 
end  surface  of  the  new  thin-slab  are  calculated  by  (31)  and 
(37),  respectively.  This  marching  algorithm  includes  all  the 
diffraction  and  forward  scattering  effects.  The  wave-type 
conversion  is  expressed  by  the  terms  Usp  and  XJPS .  Note  that 
for  this  wide-angle  forward  scattering  algorithm,  the 
wavenumber  coupling  integrals  (33)  are  not  convolutions. 
Therefore  the  corresponding  operations  in  space  domain  are  not 
local  This  shows  that  for  wide-angle  propagation,  the  correct 
scattering  effect  cannot  be  modeled  with  standard  phase- 
screens.  We  will  call  the  above  version  of  our  formulation 
represented  by  (31),  (37),  and  (38)  as  the  " thin  slab ”  method. 
The  thickness  of  the  "thin-slab"  should  satisfy  two  criteria.  One 
is  from  a  requirement  of  the  Bom  approximation.  In  order  to 
keep  the  scattered  field  small  compared  with  the  incident  field, 
roughly  speaking,  the  phase  deviation  due  to  passage  of  wave 
through  the  thin-slab  should  be  less  than  1  rad,  similar  to  equa¬ 
tion  (15)  of  Wu  [1989/?].  For  P-P  scattering,  it  is 


^-^-kaAx<\ 

oto 

where  ka  is  the  maximum  wavenumber  of  the  wave  field.  The 
other  criterion  is  that  the  propagation  within  the  thin-slab 
should  satisfy  geometric  optics,  which  requires  that  the  Fresnel 
radius  on  S j  due  to  a  propagation  distance  Ax  be  smaller  than 


the  average  size  of  heterogeneities  on  Si.  Detailed  discussion 
will  appear  in  future  publications. 

Now  we  estimate  the  amount  of  computation  involved.  The 
thin-slab  method  is  a  wavenumber  domain  marching  algorithm. 
The  model  perturbations  5p,  8A,  and  5ji  can  be  transformed  to 
the  spectral  domain  using  FFT  before  propagation.  The  compu¬ 
tation  for  this  transformation  is  negligible  compared  with  the 
propagation  calculation.  For  each  step  forward  (propagation  of 
Ar),  we  need  many  mappings  (matrix  multiplications)  from  the 
Kr  plane  to  K'r  plane.  For  Up/\  Uw,  and  U5/>,  we  need  one 
scalar  mapping  for  each.  Uss  has  two  orthogonal  components 
and  therefore  needs  two  mappings.  Therefore,  there  are  a  total 
of  5  mappings  for  each  step.  The  amount  of  computation  is 
similar  to  FM91’s  vector  phase-screen  algorithm  in  terms  of 
order  of  magnitude,  though  the  thin-slab  method  has  more  pro¬ 
jection  calculations.  However,  since  the  slab  thickness  is  only 
required  to  be  small  enough  for  the  Born  approximation  to  be 
valid,  it  can  be  larger  than  the  step  length  of  the  screen  method 
which  makes  an  additional  parabolic  approximation,  as  will  be 
discussed  in  next  section.  For  3D  modeling  the  thin-slab 
method  is  still  computationally  intensive,  although  it  is  much 
faster  than  the  3D  finite  difference  method.  In  the  following 
we  will  make  further  approximations  and  derive  a  much  faster 
elastic  complex-screen  (ECS)  method. 

4.  Small-Angle  Scattering  Approximation: 

A  Complex  Screen  Formulation 

The  forward  scattering  calculations  in  (33)  and  (31)  are 
valid  for  wide-angle  scattering.  If  we  consider  only  small-angle 
forward  propagation,  then  the  parabolic  approximation  can  be 
used  to  further  simplify  the  formulas.  Here  we  follow  the 
approximations  made  by  Wales  and  McCoy  [1983]  in  their 
derivation  using  displacement  potentials.  The  goal  is  to 
"compress"  the  thin-slabs  into  "screens",  then  replace  the 
wavenumber  domain  matrix  operations  with  space  domain 
point-to-point  multiplications,  as  in  the  scalar  wave  phase- 
screen  method.  In  this  way,  a  dual-domain  technique  using 
FFTs  can  be  applied  to  speed  up  the  computation  substantially. 

In  the  wavenumber  domain  the  parabolic  approximation  can 
be  expressed  as 

K 2 

Y=  =  *(1-  -L)  \KT\«k.  (39) 

2k 

Now  we  make  approximations  to  the  exchange  wavenumbers: 
k'a  -  ka  =  (pfl iTJr?  -  V^a  -Kf)  4  +  K'r  -  Kr 

=  (|j  •  ^jjT-)**  +  K'r  - Kr  =  04  +  K'r  -  Kr  .  (40) 

Here  we  drop  the  term  (Kjf2k2 -.K'T2l2k2),  since  it  is  generally 
a  higher  order  small  quantity  than  the  second  order  which  we 
retain.  In  a  similar  way,  we  obtain 

k'p  -  k a  ~  (&p  -  k  a)  $x  +  KV  -  K7- 

k'a  -  kp  =  ( ka  -  A'p)  ex  +  K'7-  -  Ky 

k'p  -  kp  =  0  er  +  K'r  -  K7  .  (41) 

The  3D  Fourier  transforms  of  the  perturbations  then  can  be 
approximated  by 
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dx  Sp(x,xT)e (K'r  Kt)’  %t 

x,K't  -  Kr)  =  5p(0,K'r-Kr) 

Ax 

5p(k'p  -  ka)  =  |d!xe~'<‘P*0%(x,K,J--Kr) 

=  5p(0,KV-Kr)n(A*) 

5p(k'a  -  kp)  =  6p(0,K'r  -Kr)ri*(Ar) 

8p(k'p  -  kp)  =  8p(0,K'j-  -  Kj-)  (42) 

where 

Ax 

8p(0,Kr)  =  j8p(x,Kr)<£c  =  Ax8p(Kt-)  (43) 

^(kn-  l-)Ax 

ri(Ax)  =  sinc((kfi’ka)Ax/2)e  2  p  (44) 

with  Sp(Kr)  as  the  2D  Fourier  transform  of  5p(x,xr)  with 
respect  to  xr,  and  r\*  as  the  complex  conjugate  of  T|.  Similar 
expressions  can  be  derived  for  8A.  and  8p. 

The  scattered  fields  under  this  small-angle  approximation 
become 

U"(K'r,Kr)  =  -JL-k2ua°k'a 

2ya 

8p(0,Kr)  SX(0,Kt)  28p(0,Kr) 

Po  +  2|io  Xo  +  2p<) 

=  -ikaAx  u  a°(Kr )  ------ k  'g  (45) 


UM(K'r,Kr)  =  -iLkfAx  [up(KT-)  -  £'p(up*  £'#] 

^  a 

8p(0,Kt-)  SptO.K;-) 

Po  Po 

=  - ikpAx  [up°(Kr) - k'p(u$-  -f'p)]-^2  (48) 

Po 

where  =  KV  -  is  the  horizontal  exchange  wavenumber 

and  8a,  Sp,  are  the  P  and  S  wave  velocity  perturbations, 
respectively.  We  see  that  for  the  non-coupling  P  and  S  scat¬ 
tered  waves,  the  major  contributions  come  from  the  P  and  S 
wave  velocity  perturbations.  However,  for  wave-type  conver¬ 
sion  between  P  and  S ,  only  8(3  and  Sp  have  contributions. 
The  cross-coupling  term  between  two  orthogonal  S -components 
in  (31)  has  been  dropped  in  (48),  since  it  is  a  higher  order 
small  quantity  in  forward  scattering  problems.  It  is  known  that 
the  effect  of  shear  modulus  perturbation  for  S  wave  incidence 
is  similar  to  a  double-couple  secondary  source,  which  has  its 
maximum  amplitude  of  cross-coupling  in  the  direction  perpen¬ 
dicular  to  the  incident  direction  [see  Wu  and  Aki ,  1985;  Wu, 

1 989Z?  ].  For  this  reason,  the  cross-coupling  can  be  neglected  in 
the  complex-screen  formulation  which  is  a  formulation  of 
small-angle  scattering  approximation.  In  order  to  treat  the 
cross-coupling  problem  rigorously,  we  have  to  go  back  to  the 
more  accurate  thin-slab  formulation. 

From  (45)  to  (48)  we  can  prove  that  the  scattered  waves  are 
in  the  form  of  convolution  integrals  of  the  incident  field  and  the 
perturbations  in  the  wavenumber  domain.  Therefore  in  the 
space  domain  the  perturbations  will  be  local  operators,  similar 
to  passing  "screens".  Let  us  take  P-P  and  P-S  scattering  as 
examples  to  demonstrate  how  (45)-(48)  can  be  related  to  space 
domain  screen  equations  for  the  common-mode  scattered  waves 
and  converted  waves. 

The  space  domain  total  P  wave  field  can  be  obtained  by 
Fourier  transforming  the  first  equation  of  (38).  For  the  P-P 
scattered  field,  from  (45)  we  have 


8p(k'a  -  kR)  =  | 


Uw(K'r,Kr)  *  -£-*J«®(KrXJ?«-fV(f«*  £'(s)] 

2Tp 


8p(0,Kr)  Po,8fl(0,fcr) 

T|(Ax)  -  2(— ) - - - T|(Ax) 


Po 


«o  Po 


Uw’(*i,X|r)  =  JLjrf 


x^jdKTVpp(K'T>KT) 


~  -  &  pAx  u°  (KrXf.  -  k'p(ka-  £'p)Jn(Ax) 
..  Sp(Kr)  ,2po  tx  8p(Kr) 

x  -  h  (  1 J 

Po  «o  po 


Us/>(K'r,Kr)  =  ~yki(af( Kr)-  k\)k\ 

^7  a 

8p(°,Kr)_i/A  S  Po  N  8p(0,Kr)  _  ' 

- - r)  *  (Ax )  -  2( — ) - T|*(Ax) 

Po  ao  po 


=  - ikaAx  (up(Kr)-  k>a)k'an*(Ax) 

8p(Kr)  2pp  8p(Kr) 

Po  ^  oto  Po 


(46) 


-jdK'-  VKV  f 


4k2 


re 


ikak'a 


x-^-pK T  u°(KtMK't  -  Kt)Ax 


(49) 


where  G  =  8a/a0  and  Ax  =  xj  -x.  Let 


0(x,K'r)  =  -~;jdKTu^(KT)Q.(K’T-KT)Ax  . 

Then 

VPF(X1,XIT)=  -V|d>(X|,Xi7 r) 

=  -2V|p5  <&(x  ,xT)-^-gp(x  uxlT;x  ,xT) 

=  -2V ,  [c/S  ika$a(x ,xT )&(x ,xT)Ax-^-gp (x\,xlT',x ,xT)  (50) 


x 
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where  stands  for  the y  opera1'00  with  respect  to  x,,  uw(x,ptlr)  =  —  ’{dK'Ie*7’  \dKT  Uw(KV,Kr) 

.  n  Un  47T2  J  J  * 


<F  =  “7T  »  ancl  £7  is  the  scalar  Green’s  function  with  the  P 

IK  a 

wave  propagation  speed,  defined  as 


47T 

4K2kn 


Jd  K' 


/K'. 


r  e 


T  IT  p11  p 


£P(xi;x)  = 


1 


47tjxi  -x) 

The  primary  field  can  be  written  as 


exp^ajxi-x])  . 


x[(ik^)xjdKT  ua°(Kr)S  (K'r  -  Kr)Ar 


=  prV,xV,xJ<fS  (— £^p)V4>0(jc  ,xr )B  (x  ,xr)Ax 


uS(^i.Xi7-)  = -^-J^Kre  ",rc  "a  uS(x,Kt-) 

V7a“i*,^a<|>V,Kr) 


Xgj^(^i.xi7-U.xr) 


(57) 


=  2V||j5  . 


Therefore  the  total  field 

uW*(xl>xir)  =  ua(xl  jx17*)  +  (x  i»Xj7*) 


(51) 


=  2V}\dS  <j)°(x,Xr)[l  -z£a&(x,xr) Ax] 


S 

x^pC-^hXjr^^)  . 


(52) 


where  gf  is  the  scalar  Green’s  function  with  the  S  wave  propa¬ 
gation  speed.  We  see  that  the  P-S  conversion  is  produced  by 
the  distortion  of  the  curlless  field  u£  =  V<j>°  after  passing 
through  a  screen  #(x,xr).  Since  T|(Ax)  is  a  complex  number, 
we  call  the  screen  a  "complex-screen’. 

The  calculation  would  be  inefficient  in  either  the  space 
domain  or  the  wavenumber  domain  alone.  The  ideal  way  is  to 
propagate  in  wavenumber  domain,  but  interact  with  screens  in 
space  domain,  and  shuttle  between  these  two  domains  using 
FFT.  This  is  the  dual-domain  approach.  In  the  following  we 
change  (45)  -  (48)  into  a  dual-domain  formulation. 

Dual-Domain  Implementation 
For  P-P  scattering,  from  (53)  we  have  the  total  field 

Ha** 


When  Ax  is  not  small,  it  is  better  to  use  the  Rytov  approxima¬ 
tion  instead  of  the  Bom  approximation.  After  doing  the  Rytov 
transformation  [see  Ishimaru ,  1978;  Wu  and  Flatted  1990],  we 
have 

upp(x\tXiT)  =  ~-Vj w £ (x  jXf) exp[-z£aQc(x ,Xj)Ax ] 


upp(xh K'r)  =  k'ae  a 


xT  e 


|KV"  XT 


rWa(x,xr)expH'/:aa(x,Xr)Ax]  .  (58) 


x^Sp0cuxlT;xtxT) . 


(53) 


It  can  be  seen  that  the  wave  field  interacts  with  the  screen  in 
the  space  domain,  then  is  transformed  back  to  the  wavenumber 
domain  and  propagated  to  a  new  surface  at  xi. 

For  P-S  scattering,  we  can  rewrite  (57)  as 


U«(x,,K'r)  =  (59) 


We  can  see  that  the  term  exp[-i/:aa(x,Xr)Ax]  is  the  phase- 
screen  term.  In  the  case  of  weak  distortion,  the  gradient  of  the 
wave  field  can  be  approximated  by  taking  the  gradient  of  only 
the  phase  term,  resulting  in 


where 


u£Xx,K.'T)  =  jdxTe‘KT  *Tug(x,xT)CB(x,xr)  (60) 


uw’(xi,Xir)  =  u£(x ,xr) exp[-ifca&(x ,xr)Ax ] 

> 


(54) 


is  the  distorted  P  incident  field  after  passing  through  the 
complex-screen: 

CB  (x  ,xT)  =  -iktf\(te)\dx'  B{x\xt) 

X 

where  p  is  the  unit  vector  in  the  P  wave  propagation  direction.  ifcpT|(Ax)AxZ?(x,X7)  .  (61) 

In  fact,  (54)  is  in  the  form  of  the  scalar  Rayleigh  integral  of  the  The  complex-screen  in  (61)  can  be  also  replaced  by  a  complex 

second  type,  representing  the  free  propagation  from  plane  S  to  phase-screen: 
plane  Si  of  the  distorted  wave  field.  The  distortion  is  caused 
by  a  phase-screen  with  the  amount  of  distortion  given  by  - 
kaS t(x,xr)Ax. 

For  the  P-S  conversion,  (46)  can  be  rewritten  as 


Hb (x,xt)  -  exp [—ik pr| (Ax )AxB (x ,x r ) ]  .  (62) 

This  is  because 
1 


Uw(K7,Kr)  =  /*p£'fix(fc'px£0)u°(K7-)B  (KV  -  Kr)Ax  (55) 
where 


471" 


b(kt)  =  n(Ax)[-^d.  +  (3P» .  i)M; 

Po  Oo  Ho 


(56) 


\dK’Te  T‘x'Tk\x 

k' pXe  jdxTe  T  Tu£(x,xT ) 

-V |XV iXu^U i.Xir)  =  0  . 


Therefore 


(63) 
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For  computation,  (59)  can  be  also  written  as 
U«(*„K'r)  = 

x[u<V,K'r)  -  Pp(u<l>(*,K'r)*  *'p)l  (64) 

Theoretically,  the  replacement  of  the  complex-screen  C5(x,xr) 
by  a  phase-screen  HB(x,xT)  =  exp[CB(x,xr)]  will  have  no 
effect  on  the  P-S  scattered  field.  However,  due  to  the  finite¬ 
ness  of  the  Fourier  representation  of  the  field,  the  phase-screen 
calculation  may  produce  larger  errors  than  the  complex-screen. 
In  the  following,  we  keep  the  conversion  terms  as  produced  by 
complex-screens. 

In  the  same  way,  the  S-P  scattered  field  can  be  calculated 
by 

lP0c„K'r)  =  V)-  k'a)k'a  (65) 

where 

u|i)(i,K'r)  =  |4xrc'Kr  *rup(x ,xT)C*B (x ,xT)  (66) 

with 

C*B(x,xT)  =  -ikaT\*(Ax)jdx'B(x',xT) 

X 

~  -ik  aT|  *  (Ax ) A xB  (x  ,xr ) .  (67) 

Similar  to  the  case  of  P-S,  the  complex-screen  in  (66)  can  be 
replaced  by  a  complex  phase-screen: 

H*b(x,xt)  =  exp[~/^ar|*(Ax)Ax5(x,X7*)]  (68) 

because 

-f\dK'T  e'K  r  X'Tk'a[k'a'  eiat*  \dxT  e  'i'T  ^u^x^r)] 

=  7j-V,[V1.  u?(x„x17-)]  =  0.  (69) 

For  S-S  scattering,  from  (48)  the  in-plane  S-S  scattered 
field  can  be  rewritten  as 

U“(x,,K'r)  =  -e'^r'pX^'pXuf^i.KV)]  (70) 

where 

U|P(x,K'r)  =  jdxT  e’*7  Xru^(x,xr)Cp(x,xr)  (71) 

is  the  distorted  S  incident  field  after  passing  through  the 
comp  lex -screen: 

xi 

Cp^.xj-)  =  -ikpjdx'$(x',xT)  =  -i*pArji(jE,xr)  (72) 

X 

where  [J  =  6p/p0.  Since  the  incident  S  wavefield  can  be  writ¬ 
ten  as 

UfVbK'r)  =  -e'^r'pXfr'pXupVi-KV)] 

=  -  e  (.'$x[k'$x^dxT  e  T  >'7'up(x,xr)]  .  (73) 

Then  the  total  S-S  field  can  be  written  as 


U®(*,,KV)  =  upV„KV)  +  UM(x,,KV) 

=  -e'^^ £'pX[£'pxjdxT  e'Kr  ”Tu$(x  ,xT)H  p(x  ,xT)] 

=  e  '  V*  [if >(*  ,K'r )  -  r 'p(uf  (x  ,K'r  )•  k'0]  (74) 

with 

//p(x,xr)  =  exp[Cp(x,x7-)] .  (75) 

Computation  Speed 

Let  us  estimate  the  increase  of  computation  speed  of  the 
complex-screen  method  versus  the  thin-slab  method.  From  the 
above  formulas,  we  see  that  for  one  step  forward,  the  screen 
method  needs  7  inverse  2D  FFTs,  one  scalar  for  «a(Kr),  two 
vectors  for  ua(Kr)  and  up(Kr),  and  10  forward  2D  FFTs,  one 
scalar  for  upp ,  3  vectors  for  XJPS y  U5P,  and  u55.  Each  2D  FFT 
needs  2N2l\og2N2  operations.  In  the  space  domain,  the  interac¬ 
tion  between  the  wave  field  and  the  screen  needs  only  N2 
operations.  In  the  case  of  the  thin-slab  method,  we  need  5 
mappings  (matrix  multiplications)  in  wavenumber  domain  for 
each  step  forward.  Each  mapping  needs  N2xN2  operations  if  all 
the  wavenumbers  are  taken  into  account.  Therefore  the  relative 
speed  factor  of  the  two  methods  is  about 

Sf  =  5N4/(\lN2Log2N2  +  1(W2) . 

When  N  =  128,  SF  =  330,  while  for  N  =  512,  SF  =  4148. 
For  large  3D  models,  the  complex-screen  method  can  be  2  to  3 
orders  of  magnitude  faster  than  the  thin-slab  method.  Of 
course,  the  amplitude  information  for  large  angle  scattering 
from  the  complex-screen  method  is  not  as  accurate  as  from  the 
thin-slab  method.  If  the  scattering  and  propagation  problem  is 
for  large-scale,  smoothly  varying  inhomogeneous  media,  the 
complex-screen  method  is  preferable  in  view  of  the  huge  gain 
in  computational  speed  of  the  method. 

Numerical  Examples 

To  test  our  algorithm  for  the  complex- screen  method  and  to 
compare  our  method  with  the  vector  screen  formulation  of 
FM91,  we  present  in  this  paper  two  3D  numerical  examples. 
These  examples  serve  only  to  demonstrate  the  principle  and 
capability  of  our  method.  Detailed  comparisons  with  other 
methods  and  discussion  of  accuracy  and  limitation  of  the 
method  are  deferred  for  future  publications.  Our  formulation 
has  a  major  difference  from  FM91.  In  our  theory  the  P-S  and 
S-P  conversions  are  generated  by  S  wave  velocity  perturba¬ 
tion  and  the  shear  modulus  perturbation.  For  a  Poisson  solid,  S 
wave  velocity  perturbation  plays  a  major  role  for  the  P-S  and 
S-P  conversions  under  this  parabolic  approximation.  In  the 
special  case  of  cxq  =  2(30,  S  wave  velocity  perturbation  becomes 
the  only  factor.  In  contrast,  the  P-S  conversion  of  FM91  is 
determined  by  P  wave  velocity  perturbation,  while  the  S-P 
conversion,  by  S  wave  velocity  perturbation,  which  can  cause 
nonphysical  conversions.  This  can  be  seen  in  the  following 
examples. 

I.  The  case  of  only  P  velocity  perturbation.  The  first  exam¬ 
ple  is  the  case  in  which  only  the  Lam  constant  X  is  perturbed. 
In  this  case,  there  should  be  no  P-S  conversion,  as  correctly 
predicted  by  this  theory  and  the  scattering  theory.  However, 
from  FM91,  P-S  conversion  would  still  be  generated  and  no 
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Fig.  2.  The  geometry  of  the  numerical  examples.  A  plane  P  wave  is 
incident  on  a  uniform  sphere  of  radius  R .  The  background  velocities 
are  P  velocity  Oq,  S  velocity  po-  The  parameters  in  the  two  examples 
are  (example  1)0 ($  =  4.206  km/s,  po  =  2.664  km/st  R  =7  km,  with 
10%  P  velocity  perturbation  and  (example  2)  =  6.0  km/s,  po  =  3.5 

km/s,  R  =  1  km,  with  10%  S*veIocity  perturbation. 


difference  could  be  distinguished  between  X  and  \i  perturba¬ 
tions.  Our  theory  is  supported  by  3D  finite  difference  calcula¬ 
tions.  The  experimental  geometry  is  shown  in  Figure  2.  The 
background  medium  has  P  wave  velocity  of  4.206  km/s,  S 
wave  velocity  of  2,664  km/s,  and  density  of  2.14  glcm 3.  The 
model  is  a  solid  sphere  with  10%  P  velocity  perturbation  (a 
fast  sphere)  and  no  S  velocity  or  density  perturbations.  The 
radius  of  the  sphere  is  7  km.  A  plane  P  wave  is  incident  on  the 
sphere  along  the  jt-axis  and  the  receiver  line  parallel  with  z- 
axis  is  located  8  km  from  the  center  of  the  sphere  (see  Figure 
2).  There  are  25  receivers  on  the  receiver  line  with  interval  of 
0.5  km.  The  source  time  function  is  of  Kelly  type  with  center 
frequency  at  1.0  Hz.  In  Figure  3  a  comparison  between  the 
results  of  3D  finite  difference  (FD)  and  the  elastic  complex- 
screen  (ECS  )  method  is  shown.  The  3D  FD  has  model  param¬ 
eters  of  nx  -  128,  ny  =  96,  nz  =  96,  dx  -  dy  =  dz  ~  0.25  km, 
nt  =  1500  and  dt  =  0.005  s.  The  calculation  is  performed  on 


the  128-node  nCUBE  computer  at  the  Geophysical  Center  for 
Parallel  Processing,  part  of  the  Earth  Resources  Laboratory  of 
MIT  by  C.B.  Peng.  The  CPU  time  for  this  example  is  about  29 
min.  Note  that  due  to  the  symmetry  of  the  problem,  only  one 
quarter  of  the  model  space  is  taken  into  FD  calculation.  For 
more  general  problems,  no  such  saving  of  time  can  be 
achieved.  The  complex-screen  calculation  is  performed  on  the 
SUN  SPARC  station  2  computer  at  the  Institute  of  Tectonics, 
University  of  California,  Santa  Cruz.  For  this  example,  we  use 
a  64  x  64  grid  for  the  y~z  plane  and  variable  step  length  in  x- 
direction  (0.5  km  inside  the  sphere).  The  CPU  time  is  less  than 
30  min  on  the  SPARC  station  2.  In  Figure  3,  the  distortion  of 
the  incident  P  wave  and  the  diffracted  P  waves  from  the  top 
and  bottom  points  of  the  sphere  are  seen  clearly.  There  are  no 
converted  S  waves,  consistent  with  the  scattering  theory.  It  is 
seen  also  that  the  complex-screen  method  agrees  well  with  the 
finite  difference  calculation.  In  Figure  4 a  are  shown  ECS 
snapshots  for  the  z  -component  of  the  wave  field,  which  is 
composed  of  only  scattered  waves,  where  the  generation  and 
wave  front  propagation  of  the  scattered  P  waves  are  clearly 
seen.  In  Figure  Ab  some  of  the  corresponding  snapshots  from 
FD  (A,  t  =  1  s;  B,  /  =  3  s;  C,  t  =  5  s;  D,  t  =  7  s)  are  shown 
for  comparison.  We  see  that  except  for  the  lack  of  backscat- 
tered  waves  in  the  elastic  complex-screen  method,  the  forward 
wave  fronts  are  quite  similar  for  these  two  methods. 

2.  The  case  of  only  S -velocity  perturbation.  In  this  case,  the 
major  purpose  is  to  check  the  PS  converted  waves.  Since  for 
the  FD  method  it  is  difficult  to  separate  P  and  S  waves,  we 
compare  our  results  with  an  exact  solution  for  a  solid  sphere 
with  10%  S -velocity  perturbation.  The  experimental  geometry 
is  similar  to  Figure  2,  but  with  different  parameters.  The  back¬ 
ground  medium  has  P  wave  velocity  of  6.0  km/s,  S  wave  velo¬ 
city  of  3.5  km/s,  and  density  of  2.7  glcm3.  The  radius  of  the 
sphere  is  1  km.  A  plane  P  wave  is  incident  on  the  sphere 
along  the  x-axis  and  the  receiver  line  is  parallel  with  the  z- 
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Fig.  3.  Comparison  of  3D  elastic  complex-screen  (ECS)  method  with  3D  finite  difference  (FD)  for  an  elastic  sphere  with 
10%  P  velocity  perturbation.  The  receiver  line  is  8  km  from  the  center  of  the  sphere  (see  Fig.  2).  Shown  on  the  left  are  finite 
difference  results  and  on  the  right,  ECS  results.  The  upper  part  are  synthetic  waveforms  of  the  x  -component  and  the  lower  part 
is  for  the  z  -component. 
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Fig.  Aa .  Snapshots  from  ECS  for  the  z  -component  of  scattered  P  waves  in  example  1.  The  position  of  the  sphere  is  outlined 
in  each  snapshot.  The  generation  and  propagation  of  the  scattered  waves  can  be  seen  clearly. 


axis,  located  2  km  from  the  center  of  the  sphere  (see  Figure  2).  with  frequency  interval  0.25  Hz  for  the  exact  solution,  and  0.5 

There  are  24  receivers  on  the  receiver  line  with  intervals  of  0.1  Hz  for  the  complex-screen  method.  The  exact  solution  is  cal- 

km.  The  source  time  function  is  a  Berlage  impulse  culated  using  the  eigenfunction  expansion  series  [. Petrashen , 
/(f)  =  w  exp(-w/2.5)sin  w,  where  w  =  2nf  q! ,  with  center  fre-  1945;  Ying  and  Truell ,  1956;  Korneev  and  Johnson ,  1993]  by 

quency  / o  at  10.0  Hz.  The  frequency  range  is  from  0  to  64  Hz  V.A.  Korneev  on  a  Sun-Solboume  computer  at  the  Center  for 


u  ll  if*  :  l  2ft  '  l  v,  l  ft  II  1  ft  2 1 


Fig.  Ab .  Some  corresponding  snap  shots  from  FD:  A,  t  =  1  s;  B,  t  =  3  s;  C.  t  =  5  s;  D,  t  =  7  s.  From  comparison  of  Figure 
Aa  and  Ab  it  can  be  seen  that  except  for  the  lack  of  backscattcrcd  waves  in  the  ECS  method,  the  forward  wavefronts  are  quite 
similar  for  these  two  methods. 
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Fig.  5.  Comparison  of  3D  elastic  complex-screen  (ECS)  method  with  exact  solutions  for  an  elastic  sphere  with  10%  S  velo¬ 
city  perturbation.  The  receiver  line  is  2  km  from  the  center  of  the  sphere  (see  Figure  2).  Showing  on  the  left  are  exact  solu¬ 
tions,  and  on  the  right,  ECS  results.  The  upper  part  are  synthetic  waveforms  of  the  x-component  and  the  lower  part  is  for  the 
Z  -component. 


Computational  Seismology,  Lawrence  Berkeley  Laboratory, 
University  of  California,  Berkeley.  Both  P  and  S  scattered 
waves  should  be  generated  for  a  P  wave  incidence  in  this  case. 
Here  we  compare  only  the  P-S  converted  waves.  In  Figure  5 
the  comparison  of  x-  and  z  -components  calculated  by  exact 
solution  and  by  ECS  is  given.  We  see  the  two  results  are  in 
very  good  agreement.  The  first  S  arrival  is  the  P-S  scattered 
wave  produced  by  the  right  half  of  the  spherical  surface.  The 
second  and  third  arrivals  are  the  P-S-S  scattered  waves  pro¬ 
duced  by  the  left  half  of  the  spherical  surface.  We  see  that  not 
only  the  direct  P-S  converted  waves  but  also  the  converted 
P-S  head  waves  can  be  modeled  correctly  by  the  complex- 
screen  method.  In  contrast,  there  will  be  no  scattered  waves  at 
all  from  the  formulation  of  FM91  in  this  special  case,  showing 
again  that  the  conversion  is  not  correctly  handled  in  their  for¬ 
mulation. 

The  second  difference  between  our  complex-screen  method 
and  the  vector  phase-screen  scheme  of  FM91  is  that  the 
"complex-screen"  for  converted  waves  causes  not  only  phase 
distortion  but  also  amplitude  change  of  the  incident  field.  The 
conversion  coefficient  is  Ax  dependent  as  a  sin xfx  type  func¬ 
tion.  Only  when  Ax  is  much  smaller  than  the  shear  wavelength 
does  r|(Ax)  approach  a  real  number  close  to  unity.  We  have 
proved  that  the  complex-screens  for  converted  waves  can  be 
made  into  complex  phase-screens,  but  is  not  necessary. 

For  S-S  forward  scattering  under  the  parabolic  approxima¬ 
tion  the  in-plane  scattered  waves  are  generated  by  S  wave  velo¬ 
city  perturbations,  similar  to  the  case  of  P-P  scattering  caused 
by  P  velocity  perturbations.  We  do  not  test  it  here. 

5.  Limiting  Case  of  Scalar  Waves 

In  the  case  of  scalar  waves,  the  phase-screen  formulas  have 
been  derived  through  various  approaches.  The  normal  pro*- 


cedure  of  neglecting  the  higher  order  derivatives  in  the  propa¬ 
gation  direction  leads  to  the  standard  PE  (parabolic  equation- 
type  phase-screen  [ Tappert ,  1977;  F latte'  and  Tappert ,  1975; 
Martin  and  Flatted  1988],  which  can  be  also  derived  by  the 
truncated  Taylor  expansion  of  the  square-root  operator.  Based 
on  the  symmetric  splitting  of  the  square-root  operator,  a  wide- 
angle  phase-screen  formulation  was  obtained  [Feit  and  Fleck , 
1978;  Thomson  and  Chapman ,  1983],  in  which  the  square-root 
operator  is  split  into  a  propagation  term  and  an  interaction 
term.  The  propagation  term  is  a  wide-angle  free  propagator, 
while  the  interaction  term  is  a  multiplication  by  (rc-1),  where  n 
is  the  refractive  index.  This  wide-angle  phase-screen  method 
has  better  accuracy  than  the  parabolic  approximation.  The  other 
approach  for  improving  the  accuracy  is  to  match  the  ray  equa¬ 
tion  of  the  phase-screen  method  which  is  a  one-way  equation, 
with  the  ray  equation  of  the  two-way  Helmholtz  equation  [Tol¬ 
stoy  et  al.9  1985;  Berman  et  al,  1989],  In  Tolstoy  et  al  [1985], 
only  the  interaction  term  is  changed  into  lo gn,  while  the  pro¬ 
pagation  term  remains  as  the  standard  parabolic-type  propaga¬ 
tor;  in  Berman  et  al.  [1989]  both  the  interaction  and  propaga¬ 
tion  terms  are  changed.  In  this  section,  as  a  limiting  case  of  our 
wide-angle  one-way  elastic  wave  equations,  we  derive  a  thin- 
slab  equation  for  the  scalar  case,  which  may  be  termed  as  a 
generalized  wide-angle  phase-screen  equation.  The  equation 
matches  the  scattering  in  the  forward  direction  with  the  two- 
way  wave  equation.  The  method  should  be  more  generally 
valid  than  the  other  phase-screen  methods,  since  it  only 
requires  the  weak  scattering  inside  the  thin  slab  of  thickness 
Ax,  and  makes  fewer  approximations  than  the  other  methods. 
After  taking  the  small-angle  approximation  for  the  interaction 
term,  the  thin-slab  equation  degenerates  to  the  wide-angle 
phase-screen  equation  of  Feit  and  Fleck  [1978]. 

From  equation  (30)  and  (31),  consider  only  the  P-P  scatter¬ 
ing,  and  set  |i0  =  0,  5p  =  5p  =  0.  The  scattered  displacement 
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field  becomes 


U(KY,Kr)  =  (r-Y.KV  -KT) .  (76) 

Y  ccq 


<K*i>Xir)  =  2^dS  |<|»0(x ,Xt ) - ik <t>°(x ,xT )0t(x ,xr)A*j 


The  total  field  is 


x—g(xuxlT-x,xT) 


(82) 


u(x,,K'r)  =  e^xk' 


or 


i 


u°(x,K't)  -  iyfdKr  «°(Kr)-^-(Y'  -  Y>  K'r  -  Kr) 


xir)  =  2|<fS  <J)°(;c,X7-)exp 


-ikO.(x  ,xT)  Ax 


•  (77) 


x-^g(xi,xiT;x,xT) . 


(83) 


Set 


u(Kr)  =  ikic$(KT) 


where  (J>(Kr)  is  the  displacement  potential  which  satisfies  the 
scalar  wave  equation,  we  have 


<K*I,KV)  =  e1^ 


Ux ,K'r)  - i^-U Kt  <j>0(x ,Kr)— (/ - y,  K'r  - Kr) 

Y  J  Oo 


1  -i—r  IdK 

i J 


=  e^UxXr) 

,Kr)  §g 
4>o(x  ,K'r)  c<o 


(Y'-Y.  KY  -Kr) 


(78) 


Assuming  that  Ax  is  made  up  of  many  tiny  steps,  the  above 
expression  can  be  put  into  a  form  of  phase-screen: 


<!>(*,, KY)  =  e** 

where 

A\)/(Ax,K'r)  =  -~J<7K 


<!>o(x.K'r)e 

(j>o(x,Kr)  §Qf 
<})o(x,KV)  oto 


-/Ay(Ax,K'  } 


(79) 


(Y  -y,  K't  -Kr)  (80) 


Equations  (78)  and  (79)  are  the  wide-angle  thin-slab  formulas 
in  the  wavenumber  domain.  For  the  wave  field  at  each 
wavenumber,  the  phase  should  be  distorted  before  propagation. 
The  amount  of  distortion  depends  on  the  spectrum  of  the 
heterogeneities  within  the  slab.  However,  the  distortion  is  not  a 
local  operation  in  the  wavenumber  domain.  Instead,  the  amount 
of  distortion  of  the  field  at  one  wavenumber  is  the  summation 
of  the  contributions  from  the  incident  field  at  all  the 
wavenumbers  interacting  with  the  heterogeneities.  Since  the 
interaction  between  the  incident  field  and  the  heterogeneities  is 
not  a  convolution  in  the  wavenumber  domain,  the  operation  in 
space  domain  is  not  local,  i.e.,  not  a  multiplication.  Therefore, 
strictly  speaking,  (78)  and  (79)  are  not  the  standard  formulas  of 
phase-screen,  but  may  be  termed  as  a  generalized  phase-screen 
method.  In  the  following  we  will  show  that  under  the  narrow- 
angle  approximation,  these  formulas  approach  the  regular  for¬ 
mulas  of  scalar  phase-screen. 

Following  a  similar  procedure  as  for  the  case  of  P-P  elastic 
wave  scattering,  (78)  can  be  simplified  using  the  parabolic 
approximation  as 


\)  =  e 


»YAx 


<t>o(x  >K'r)  -  ikjd Kj  <t>o(jt ,Kr)a(x  ,K'r  -  KT)Ax 


(81) 


where  Cc(jc  ,Ky-)  is  the  2D  Fourier  transform  of  the  slab  at  x. 
Upon  transforming  into  the  space  domain  we  obtain  a  regular 
form  of  scalar  phase-screen: 


Note  that  in  the  above  derivation,  the  parabolic  approximation 
or  narrow-angle  approximation  is  only  applied  to  the  interaction 
between  the  incident  field  and  the  heterogeneities,  the  free 
space  propagation  remains  the  wide-angle  one  or  can  be 
approximated  by  any  asymptotic  form.  It  seems  arguable  that 
these  two  approximations  should  be  kept  in  the  same  order; 
however,  these  two  are  decoupled  in  some  degree.  The  interac¬ 
tion,  namely,  the  scattering  process,  strongly  depends  on  the 
properties  of  the  heterogeneities,  mainly  the  smoothness  and  the 
average  size  of  the  heterogeneities.  When  where  a  is 

the  average  size  of  the  heterogeneities,  large-angle  scattered 
waves  have  very  little  energy,  and  therefore  the  small-angle 
approximation  can  be  applied.  Of  course,  the  scattered  ampli¬ 
tudes  of  wide-angle  scattered  waves  will  carry  larger  errors 
than  the  small-angle  ones.  On  the  other  hand,  the  approxima¬ 
tion  to  the  free  propagator  depends  on  the  observation  angles. 
The  use  of  a  wide-angle  propagator  will  keep  the  correct  phase 
(or  travel  time)  information  even  though  the  amplitude  has  rela¬ 
tively  large  errors  due  to  the  small-angle  approximation  applied 
to  the  interaction.  In  fact,  (83)  is  nearly  equivalent  to  the 
wide-angle  phase-screen  formula  of  Feit  and  Fleck  [Feit  and 
Fleck ,  1978;  Thomson  and  Chapman ,  1983;  Wu  and  Huang , 
1992],  which  has  been  shown  to  have  better  performance  than 
the  narrow-angle  phase-screen  method.  However,  if  we  want 
to  adopt  both  wide-angle  interaction  and  wide-angle  propaga¬ 
tion,  the  thin-slab  method  (78)  or  (79)  has  to  be  used. 

6.  Conclusion 

Wide-angle  one-way  elastic  wave  propagation  in  arbitrarily 
heterogeneous  media  is  formulated  using  the  elastic  Rayleigh 
integral  and  elastic  Bom  scattering  theory  in  both  space  domain 
and  wavenumber  domain.  The  wavenumber  domain  formulation 
leads  to  a  compact  solution  to  the  one-way  propagation  and 
scattering  problems.  It  is  shown  that  for  wide-angle  scattering 
the  scattering  effects  of  a  thin-slab  cannot  be  equated  to  pas¬ 
sage  through  a  regular  phase-screen,  since  the  interaction 
between  the  incident  wave  field  and  the  slab  is  not  local  in  both 
the  space  domain  and  the  wavenumber  domain.  We  call  this 
more  generally  valid  formulation  the  "thin-slab"  formulation. 
After  applying  the  small-angle  approximation,  the  thin-slab 
effect  degenerates  to  that  of  an  elastic  complex-screen. 

Compared  with  the  scalar  phase-screen,  the  elastic  phase- 
screen  has  the  following  features:  (1)  For  P-P  scattering  and 
S-S  in-plane  scattering,  the  elastic  complex-screen  acts  as  two 
separate  scalar  phase-screens  for  P  and  S  waves  respectively. 
The  phase  distortions  are  determined  by  the  integrated  P  and  S 
wave  velocity  perturbations  respectively.  (2)  For  PS  and  S-P 
conversions,  the  screen  is  no  longer  a  pure  phase-screen  and 
becomes  complex  (with  both  phase  and  amplitude  terms);  the 
magnitude  and  phase  of  the  complex  factor  depend  on  the 
thickness  of  the  thin-slab,  frequency  and  the  difference  between 
the  P  and  S  wave  velocities.  Both  conversions  are  functions  of 
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the  shear  wave  velocity  perturbation  and  the  shear  modulus  per¬ 
turbation.  For  Poisson  solids  the  S  wave  velocity  perturbation 
plays  a  major  role.  In  the  special  case  of  cxq  =  2p0,  the  S  wave 
velocity  perturbation  becomes  the  only  factor  for  both  conver¬ 
sions.  (3)  For  the  cross-coupling  between  in-plane  S  wave  and 
off-plane  S  wave,  only  the  shear  modulus  perturbation  5ji  has 
influence.  For  the  complex-screen  method  the  cross-coupling 
term  is  neglected  because  the  term  is  a  higher  order  small 
quantity  for  small-angle  scattering. 

Comparisons  with  solutions  from  the  3D  finite  difference 
method  and  exact  eigenfunction  expansions  are  made  for  two 
special  cases.  One  is  for  a  solid  sphere  with  only  P  velocity 
perturbation;  the  other  is  with  only  S  velocity  perturbation.  The 
elastic  complex-screen  method  agrees  generally  well  with  the 
3D  FD  method  and  the  exact  solutions.  Comparison  of  the 
results  in  this  paper  with  that  of  FM91  showed  that  the  latter 
method  does  not  handle  wave  conversions  properly,  especially 
for  the  P-S  conversion  and  the  cross-coupling  between 
differently  polarized  S  waves. 

In  the  limiting  case  of  scalar  waves,  the  derivation  in  this 
paper  leads  to  a  more  generally  valid  new  method:  a  scalar 
thin-slab  method.  When  making  the  small-angle  approximation 
to  the  interaction  term  while  keeping  the  propagation  term 
unchanged,  the  thin-slab  method  approaches  the  currently  avail¬ 
able  scalar  wide-angle  phase-screen  method. 


Appendix  A:  Elastic  Wave  Green’s  Tensors 
in  Spase  Domain 

For  an  infinite  homogeneous  and  isotropic  elastic  medium, 
the  Green’s  displacement  tensor  (dyadic)  is 


where 


e,e/=Z«.«i  =  1  • 
1=1 


Note  that  V'G(x,x')  =  -VG(x,x'). 
Green’s  strain  tensor  (triadic): 
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Green’s  traction  tensor: 
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where  ka  =  {Q/a0,  kp  =  co/p0, _ r_=  Ix-xl,  r  =  (x-x')/r, 

«Xo  =  V(Ao  +  2pt{,)/p0,  and  (30  =  Vph/Po  are  the  P  and  S  wave 
propagation  speeds  in  the  background  medium,  respectively, 
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The  gradient  of  Green’s  displacement  tensor  (triadic) 
V'G(x,x')  =  V'Ga  +  V'Gp 
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Appendix  B:  Elastic  Wave  Green’s  Tensors 
in  Wavenumber  Domain 

In  homogeneous  media,  Green’s  displacement  tensor  can  be 
expressed  as 

G(x,x')  =  -L-  j^-p  I  +  VV)gs  (x,x')  -  VVgp  (x,x')j  (B 1 ) 
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with  ka  =  o)/a  and  k§  =  co/(3. 

In  the  wavenumber  domain,  we  know  from  the  Weyl  integral 
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From  (B3)  we  derive  the  gradients  of  gp  and  gs : 


(B4) 


V*'=8 


8tt 

then  Green’s  displacement  tensor, 


(B5) 


;*2  r 

Ga(x,xO  =  — jiyJrfKrfjfa 


87t2p(0 


*  -VV 


Ya 


ilrZ 


1  'Vr 


87C2pC0 


Yp 


its  gradients, 


(B6) 


V'Ga(x,x')  =  •  —  2  pKr^a^a — e 

87t2pco2  J 


1  'V r 


Ya 


V'Gp(x,x')  =  -Ad  Kr  £P(I  -  JfpiCp)— e,kP‘ r 

8k2pco2J  Yp 


87t2pc0' 


Ya 


V-  Gp(x,x')  =  0 
and  Green’s  stress  tensor. 


(B7) 


£(x,x')  =  XI(V-  G)  +  pi[VG  +  (VG)2 


S“(x,x')  =  -  —2^.\dKT 
Sn  pa  J 


Leik“r 

'J  Ya 


£p(x,x') 


k$ 


87T2pC0‘ 


yjika  +  2jx&akaka 

j* d  Kt  p  £pl  +  6j  £p Sj  +  2£p£p£ 


1  /yr 


X — e 
Yp 


(B8) 


Green’s  traction  tensor  on  the  surface  perpendicular  to  the  x 
axis  with  normal  along  the  positive  x  direction  (  ex  -  -n  )  is 
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A  Comparison  between  Phase  Screen,  Finite  Difference,  and 
Eigenfunction  Expansion  Calculations  for  Scalar  Waves 
in  Inhomogeneous  Media 

by  Yin-Bin  Liu  and  Ru-Shan  Wu 

Abstract  Phase  screen,  fourth-order  finite  difference  (FD),  and  eigenfunction 
expansion  calculations  of  scalar  wave  propagation  in  two-dimensional  (2D)  in¬ 
homogeneous  media  are  compared  to  assess  the  accuracy  of  the  phase  screen 
method.  The  phase  screen  method  is  a  forward  propagation  (one-way  wave) 
algorithm.  "The  finite  difference  and  eigenfunction  expansion  calculations,  which 
are  solutions  of  full  wave  equation,  are  chosen  as  references  in  this  study.  Com¬ 
parison  of  synthetic  seismograms  by  phase  screen  and  finite  difference  methods 
is  made  for  four  kinds  of  models:  (1)  multi-uniform-cylinder  model,  (2)  Gaus¬ 
sian  random  media,  (3)  exponential  random  media,  and  (4)  flicker-noise  random 
media.  Results  show  good  agreement  for  weak  random  media  (velocity  pertur¬ 
bations  s=10%).  For  discrete  heterogeneities,  such  as  the  multi-uniform-cylinder 
model,  the  results  agree  well  for  up  to  50%  deviation  in  velocities.  The  com¬ 
puter  CPU  time  of  the  phase  screen  program  for  a  problem  of  grid  size  1024  by 
512  is  367  sec  in  a  SUN  SPARC  station  II,  about  57  times  faster  than  the  FD 
program  we  used.  For  large  3D  problems  the  time  saved  is  expected  to  be  much 
greater.  For  a  single  cylinder  scatterer  with  and  without  absorption,  we  compare 
synthetic  seismograms  by  the  phase  screen  method  and  by  the  eigenfunction 
expansion  method  (exact  solution).  The  agreement  between  the  two  methods 
demonstrates  that  the  phase  screen  method  can  also  give  good  results  for  in¬ 
homogeneous  absorbing  media. 


Introduction 

Wave  propagation  in  inhomogeneous  media  has  been 
extensively  studied  (Chernov,  1960;  Tatarskii,  1961; 
Ishimaru,  1978;  Wu  and  Aki,  1985,  1988,  1989,  1990). 
In  theoretical  studies,  weak  scattering  approximations  are 
assumed  in  many  cases  to  make  the  problems  tractable. 
For  wave  propagation  in  highly  heterogeneous  media, 
however,  weak  scattering  approximations  may  not  be 
applicable.  Various  numerical  techniques,  such  as  finite 
difference,  finite  element  methods,  and  phase  screen  and 
other  one-way  propagation  methods,  are  used  in  this  case 
to  produce  synthetic  seismograms. 

The  finite  difference  method  can  produce  a  full  so¬ 
lution  of  the  wave  equation  (including  converted,  dif¬ 
fracted,  and  normal  modes).  Unlike  various  high-fre¬ 
quency  approximations  (WKBJ,  ray  method),  there  is,  in 
principle,  no  limitation  on  the  ratio  of  scatterer  size  to 
wavelength,  if  the  discretization  is  chosen  appropriately. 
The  basic  restriction  to  this  method  is  the  speed  and 
memory  size  of  the  computer.  Given  the  grid  size  con¬ 
straint  imposed  by  the  computer  speed  and  memory,  the 


accuracy  and  stability  considerations  limit  the  number  of 
wavelengths  contained  in  the  grid.  So  the  CPU  time  and 
memory  size  consideration  often  prohibit  the  use  of  fi¬ 
nite  difference  technique  for  many  outstanding  problems 
in  seismology. 

The  method  of  “phase  screen”  is  used  as  a  one-way 
propagation  method  in  problems  involving  wave  propa¬ 
gation  in  smoothly  inhomogeneous  media.  This  method 
has  been  extensively  used  to  study  light  transmission 
through  the  atmosphere  (Ratcliffe,  1956;  Mercier,  1962; 
Martin  and  Flatte,  1988),  light  signals  in  optical  fibers  (Feit 
and  Fleck,  1978),  radio  signals  through  the  ionosphere 
(Buckley,  1975;  Bramley,  1977;  Knepp,  1983),  acoustic 
waves  in  the  ocean  (Flatte  et  al.,  1979;  Thomson  and 
Chapman,  1983;  Thomson,  1990),  and  seismic  waves  in 
the  earth  (Stoffa  et  al.,  1990;  Wu  and  Huang,  1992).  Re¬ 
cently,  there  are  efforts  to  extend  the  phase  screen  method 
to  elastic  wave  propagation  (Fisk  and  McCartor,  1991, 
1993;  Fisk  et  al.,  1992;  Wu,  1994).  This  one-way  prop¬ 
agation  method  neglects  all  the  backscattered  waves,  which 
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should  be  small  compared  with  the  forward  scattered  waves 
for  large  scatters  or  smoothly  heterogeneous  media,  but 
retains  all  the  multiple  forward  scattered  waves.  There¬ 
fore,  the  method  can  model  all  the  forward  multiple  scat¬ 
tering  phenomena,  including  focusing  and  defocusing, 
diffraction  and  interference,  multi-pathing,  etc.,  with  a 
tremendous  saving  of  computation  time  and  data  storage. 

Calculating  the  propagation  and  scattering  of  acous¬ 
tic  and  elastic  waves  by  irregularly  shaped  objects  in  an 
exact  manner  is  either  impossible  or  extremely  difficult 
(Malischewsky,  1987).  However,  for  some  simple  shapes 
of  scatterers,  such  as  cylinders  and  spheres,  exact  so¬ 
lutions  have  been  given  by  solving  the  wave  equation 
and  matching  the  boundary  conditions  (Pao  and  Mow, 
1973;  Stanton,  1988).  In  order  to  assess  the  accuracy  of 
the  phase  screen  method  for  absorbing  media,  we  choose 
the  acoustic  wave  scattering  by  an  absorbing  circular  fluid 
cylinder  as  the  reference  solution. 

Numerical  Methods  and  Algorithms 

Phase  Screen  Method 

For  an  isotropic  acoustic  inhomogeneous  media  with 
propagation  velocity  v(x),  the  displacement  field  u(x) 
satisfies  the  wave  equation  (e.g.,  Wu  and  Huang,  1992) 

V2U  =  r—  u,  (1) 

v  (x) 


with 


v(x)  =  v0  +  <5v(x), 

where  v0  is  the  background  velocity,  <5v(x)  represents  the 
velocity  perturbation  relative  to  v0,  V2  is  the  Laplacian, 
and  a)  is  the  angular  frequency  of  the  wave  field.  A  time 
factor  of  e~iwt  and  a  2D  Cartesian  geometry  are  assumed. 

In  the  phase  screen  method,  the  effect  of  the  velocity 
perturbation  (absorbing  media  corresponding  to  complex 
velocity)  within  a  thin  slab  of  thickness  A z  is  treated  by 
multiplying  the  incident  wave  field  by  position-depen¬ 
dent  (in  the  2D  case,  jc-dependent)  phase  factors  elA  where 


(O 

A(x,  Az)  =  - 

V0 


(2) 


FFT  to  wavenumber  domain  and  propagated  to  the  next 
screen  in  the  z  direction;  (3)  the  displacement  field  is 
transformed  by  IFFT  back  to  space  domain  and  multi¬ 
plied  by  the  screen  factor;  (4)  the  procedures  of  (2)  and 
(3)  are  repeated  for  all  screens  until  the  observation  line 
is  reached. 

Criteria  for  Phase  Screen  Applications .  The  applica¬ 
tion  of  the  phase  screen  method  generally  requires  that 
the  screen  interval  satisfies  the  following  criteria  (Wu, 
1988):  first,  the  weak  scattering  approximations  must  be 
satisfied  for  each  screen.  Second,  the  scattering  is  pre¬ 
dominantly  of  small-angle  scattering  in  the  forward  di¬ 
rection  which  requires  ka  >  1 ,  where  a  is  the  scale  length 
of  heterogeneities,  k  =  a)/v0  is  the  wavenumber.  Finally, 
the  geometrical  optics  treatment  of  propagation  between 
screens  is  valid,  i.e.,  A  =  A z/ka2  <  1,  where  A  is  the 
diffraction  parameter  and  Az  is  the  screen  interval. 

The  computer  program  used  in  this  study  is  modified 
from  Wu  and  Huang  (1992)  to  include  the  case  of  ab¬ 
sorbing  heterogeneous  media. 

Finite  Difference  Method 

The  finite  difference  method  solves  equation  (1)  nu¬ 
merically  by  replacing  the  derivatives  in  space  and  time 
by  finite  difference  approximations.  The  real  velocity  v 
is  a  function  of  position.  At  present  there  is  no  efficient 
scheme  for  using  the  finite  difference  method  to  solve 
wave  equations  in  heterogeneous  absorbing  media. 

It  is  well  known  that  the  finite  difference  method  has 
numerical  grid  dispersion  and  stability  problems.  These 
problems  are  caused  by  the  discrepancy  between  the  fi¬ 
nite  difference  and  analytical  representations  of  both  spatial 
and  time  derivatives.  For  a  given  frequency  and  model 
velocity  distribution,  this  discrepancy  can  be  reduced  by 
decreasing  the  sampling  intervals  of  both  space  and  time. 
However,  the  reduction  of  sampling  intervals  will  dras¬ 
tically  increase  the  number  of  mesh  points  in  represent¬ 
ing  the  wave  field  in  the  spatial  and  time  domain  and 
thus  requires  substantially  more  core  memory  and  com¬ 
putation  time.  A  reasonable  compromise  is  that  at  least 
five  nodes  per  shortest  wavelength  are  needed  for  a  fourth- 
order  finite  difference  code.  The  time  sampling  interval 
can  then  be  obtained  by  the  usual  spectral  analysis.  The 
stability  criterion  can  be  stated  as  (Alford  et  al. ,  1974) 


Here,  we  assume  the  wave  is  propagating  in  the  z  di¬ 
rection  (for  derivation  of  the  phase  screen  method  see 
Thomson  and  Chapman,  1983;  Wu  and  Huang,  1992; 
Wu,  1994). 

The  phase  screen  algorithm  for  scalar  waves  may  be 
summarized  as  follows  for  the  2D  case:  (1)  given  a  2D 
grid  of  velocity  distribution,  the  medium  can  be  equated 
to  a  series  of  phase  screens;  (2)  the  displacement  field 
in  the  space  domain  along  the  x  axis  is  transformed  by 


^max  dt 
dx 


(3) 


where  vmax  is  the  maximum  velocity  in  the  model  and  dx 
and  dt  are  the  spatial  and  temporal  sampling  intervals, 
respectively. 

Boundary  Conditions .  Boundary  conditions  have  to  be 
applied  to  the  four  edges  of  the  model.  Three  types  of 
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boundary  conditions  exist.  The  first  is  Dirichlet  or  rigid 
surface  condition.  It  requires  that  the  displacement  wave 
field  at  boundary  by  zero,  i.e.,  ux  =  u2  =  0.  The  second 
type  is  the  stress-free  condition,  which  may  sometimes 
be  called  the  Neumann  or  free-surface  condition.  The 
third  type  is  the  radiation  condition,  which  may  also  be 
called  the  absorbing  boundary  or  transparent  condition. 
It  makes  the  boundary  transparent  for  incident  waves.  In 
this  study  we  choose  the  absorbing  boundary  condition 
given  by  Clayton  and  Engquist  (1977)  for  the  four  edges. 
The  absorbing  boundary  condition  can  perfectly  absorb 
the  incident  wave  perpendicular  to  the  boundary,  but  only 


(a)  Cylinder  Model 
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(b)  Multi-Uniform-Cylinder  Model 
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partially  absorb  obliquely  incident  waves.  In  order  to  re¬ 
duce  the  boundary  effects,  we  enlarge  the  model  size  so 
that  reflected  waves  from  the  artificial  boundaries  come 
much  later  than  the  waves  of  interest. 

The  finite  difference  code  used  in  this  study  is  fourth 
order  in  space  and  second  order  in  time,  developed  by 
Dr.  X.  B.  Xie  based  on  Alford  et  aL  (1974). 

Exact  Solution  of  a  Circular  Cylinder 

The  scattering  of  acoustic  and  elastic  waves  by  a 
cylinder  inclusion  has  drawn  the  attention  of  a  number 
of  authors  in  connection  with  acoustics  and  seismology 


(c)  Gaussian  Random  Media 


(e)  Flickcr-Noise  Random  Media 


Transverse  Distance 


Figure  1 .  Velocity  distributions  of  the  five  types  of  heterogeneities  considered 
in  this  study.  The  size  of  the  model  space  for  the  absorbing  cylinder  model  is 
51.2  by  51.2  km2  and  for  the  other  four  cases  is  102.4  by  51.2  km2.  Models  c, 
d,  and  e  are  constructed  from  the  same  random  number  seed.  Areas  with  higher 
than  average  velocity  arc  bright,  (a)  Absorbing  cylinder  model;  (b)  multi-uni¬ 
form-cylinder  model;  (c)  Gaussian  random  medium;  (d)  exponential  random  me¬ 
dium;  (e)  flicker-noise  random  medium. 
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Table  1 
List  of  Models 


Scale 

Correlation 

2D  Power 

Length  a 

Taper  Length 

Model 

Medium 

Function 

Spectrum 

-  \V&/  nro  (km) 

kp 

Filter 

(Grid  Points) 

la 

absorbing 

deterministic 

-0.2 

5 

10.5 

Y 

16 

lb 

cylinder 

model 

-0.2 

5 

10.5 

Y 

16 

lc 

models 

-0.3 

5 

10.5 

Y 

16 

Id 

+0.3 

5 

10.5 

Y 

16 

le 

-0.5 

5 

10.5 

Y 

16 

If 

+0.5 

5 

10.5 

Y 

16 

2a 

multi-uniform- 

deterministic 

+0.3 

4;  5 

8.4;  10.5 

Y 

16 

2b 

cylinder 

model 

-0.3 

4;  5 

8.4;  10.5 

Y 

16 

2c 

models 

+0.5 

4;  5 

8.4;  10.5 

Y 

16 

2d 

-0.5 

4;  5 

8.4;  10.5 

Y 

16 

3a 

3b 
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Note:  ko  : 

=  2tt/o/v0,  where  v0 

=  6  km/sec  is  the  background  velocity  and  /0  = 

2  Hz  is  the  central  frequency. 

Absorbing  B.C. 
t  v  -t  y  j-auuur  y 

y  Receiver 

General  Solution  of 

a  Circular  Cylinder. 

The  pressure 

X-T-3C-X-X 

field  D;  outside  (i  — 

1)  and  inside  (/  = 

2)  the  cylinder 

satisfy  the  wave  equation 


■2„  - 


V2P, 


v?  dt2  ’ 


(4) 


tilt 

Incident  Plane  Wave 

Figure  2.  Schematic  diagram  of  the  2D  grid. 
Receivers  are  located  at  the  grid  points  on  the  sur¬ 
face.  The  absorbing  boundary  condition  is  im¬ 
posed  to  the  four  edges.  Plane  acoustic  waves  are 
incident  upon  the  bottom  of  the  grid  and  propa¬ 
gating  toward  the  surface. 


where  v,(/  =  1,  2)  are  the  complex  compressional  ve¬ 
locities  outside  and  inside  the  cylinder,  respectively.  The 
boundary  conditions  require  that  pressure  and  the  radial 
component  of  the  particle  displacement  be  continuous  at 
the  boundary  of  the  cylinder. 


f?nc(a)  +  p™'(a)  =  pinl(a) 


(Pao  and  Mow,  1973;  Stanton,  1988).  The  exact  solu¬ 
tions  are  known  for  the  following  cases:  rigid,  soft,  fluid, 
and  elastic  cylinders.  In  this  study  we  consider  only  the 
scattering  of  a  plane  acoustic  wave  by  an  absorbing  fluid 
cylinder  as  a  reference  solution  to  compare  with  the  phase 
screen  method.  The  formulation  in  this  article  follows 
the  approach  of  Stanton  (1988). 


urnc(a)  +  uT\a)  =  uT(a), 


(5) 


Jnt 


where  a  is  the  radius  of  the  cylinder,  /?mc,  p" 
u'ff,  uT\  and  iff1  are  the  incident,  scattered,  and  internal 
(r  <  a)  pressures  and  radial  components  of  the  particle 
displacements,  respectively.  The  pressure  and  the  radial 
component  of  the  particle  displacement  are  related  by 
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pco2  dr 


Assume  a.  normally  incident  plane  acoustic  wave  with 
source  spectrum  g(o))  propagating  along  the  z  direction, 


g((o)ei^z~0,t)  d(o . 


The  incident  field  can  be  expanded  in  terms  of  Bessel 
function 


axis  of  the  cylinder  (0  =  0  is  the  forward  direction),  (o 
(6)  is  the  angular  frequency,  Jm(k\r)  is  the  Bessel  function 
of  the  first  kind  of  order  m.  The  general  solution  of  the 
internal  and  scattered  pressures  are 


g((o)  ^  Am{k,  w)  cos  C mO)JJk2r)e  ,u>t  do) 


g(to)  X  *0  cos  (mO)rfm\kir)e~io>t  do). 


00  00 

g((0)  ejm  cos  (mO)Jm(klr)e~‘al  dco,  (8) 

~oo  m= 0 

where  em  is  Neumann  factor  (c0  —  1  ;  2,  m  —  1,2, 

=  o)/vl  is  the  wavenumber  outside  the  cylinder, 
6  is  the  azimuthal  angle  in  the  plane  perpendicular  to  the 


where  H^Xkr)  is  the  first  kind  Hankel  function,  defined 
as  H^Xkr)  =  Jm{kr)  +  iNJJcr),  where  NJkr)  is  the  Neu¬ 
mann  function.  The  coefficients  AJJc ,  cu)  and  BJk,  o)) 
are  determined  by  boundary  conditions. 

Substituting  equations  (8)  and  (9)  into  equation  (6), 
we  can  derive  the  radial  particle  displacements, 


(a)  Multi-Uniform-Cylinder  Model  (e  =  -0.3,  a  =  5.0  km) 


(c)  Multi-Uniform-Cylindcr  Model  (e=  -0.3,  4,  16 ) 


(b)  Flicker-Noise  Random  Media  (  £  =  0.1  ,  a  =  5.0  km  ) 


(d)  Flicker-Noise  Random  Media  (£-  0.1,  a  -  5.0  km,  4,  16  ) 


Figure  3  Synthetic  phase  screen  seismograms  without  wavenumber  filter:  (a) 
for  model  2b  and  (b)  for  model  5b.  Synthetic  phase  screen  seismograms  with 
wavenumber  filter  for  different  screen  intervals:  (c)  for  model  2b  and  (d)  tor 
model  5b. 
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(a)  Multi-Cylinder  Model  (FD) 


(b)  Multi-Cylinder  Model  (Phase  Screen) 
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Figure  4.  Comparison  of  gray-scale  graphics  of  seismograms  from  finite  dif¬ 
ference  (FD)  and  phase  screen  calculations:  (a)  and  (b)  are  for  model  2a  (multi- 
cylinder  model);  (c)  and  (d)  arc  for  model  5b  (flicker-noise  random  medium). 


(10) 
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where  p,-  is  the  density  outside  (/  =  1)  and  inside  (/  = 
2)  the  cylinder  and  the  prime  on  each  J  and  H  represents 
the  derivative  with  respect  to  argument  r. 

Substituting  equations  (8),  (9),  and  (10)  into  equa¬ 
tion  (5),  we  can  solve  the  coefficients  Am(k,  w)  of  the 
internal  field  and  BJk,  w)  of  the  scattered  field 


Am(k,  b>)  = 


iraAn 


Bm(k,  co)  = 


(ID 


(c)  Flicker-Noise  Random  Media  (FD) 


Distance 


(d)  Flicker-Noise  Random  Media  (Phase  Screen) 


Distance 


Figure  4. — Continued 
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where- -  - 

A* i  =  gkxJJJc2a)Jm+x{k{a)  - 


--{g-\)Jm{kxa)Jm{k2a) 

a 

Am  =  k^m+l(k^\kxa)  -  gkJJM0&x  ( kxa ) 


+  -(g-  1) Jm(k2a)H™  ( kta ) 
a 


g  = 


P2V2 

Pifi 


(12) 


The  total  pressure  poxA  outside  the  cylinder  can  be  writ¬ 
ten  as 


p'° *  = 


g(o>) 


+  X  Bm(k,  co)  cos 

m  —  0 


•e-^dco. 


(13) 


Equation  (13)  is  the  analytic  solution  of  the  acoustic  wave 
scattering  by  a  circular  cylinder.  For  an  absorbing  cir¬ 
cular  cylinder  the  corresponding  velocity  and  wavenum¬ 
ber  become  complex.  The  incident  field  can  be  calcu¬ 
lated  by  equation  (7).  The  rate  of  convergences  of  the 
series  in  equation  (13)  depends  on  the  difference  be¬ 
tween  the  inner  and  outer  refraction  indexes,  but  is  much 
more  sensitive  to  the  value  of  kxa.  The  higher  the  values 
of  kxa  is,  the  slower  the  convergence.  In  this  article,  the 
number  of  terms  retained  in  the  series  is  60.  This  number 
is  more  than  adequate  to  render  the  truncation  errors 
negligible  in  all  the  cases  treated  in  the  article.  For  each 
frequency  loop,  we  stop  the  series  summation  when  the 
relative  difference  between  successive  terms  in  less  than 
10“5.  We  can  also  calculate  the  internal  pressure  p'n{  by 
equation  (9). 

Generation  of  Inhomogeneous  Media 

with  Absorption 

A  number  of  inhomogeneous  media  were  used  to 
compare  the  phase  screen  and  finite  difference  methods. 
We  chose  five  kinds  of  models:  (1)  absorbing  cylinder 
model,  (2)  multi-uniform-cylinder  model,  (3)  Gaussian 


Distance  (km) 


(a)  Multi-Uniform-Cylinder  Model  (  e  =  0.3,  a=  4.0,  5.0  km) 


(c)  Multi-Uniform-Cylinder  Modet  (  e  =  0.5,  a  ■  4.0,  5.0  km) 


(b)  Multi-Uniform-Cylinder  Model  (  £  =  -  0.3,  a  =  4.0,  5.0  km)  (d)  Multi-Uniform-Cylinder  Model  (  e  =  -  0.5,  a  =  4.0.  5.0  km) 

Figure  5.  Comparison  of  synthetic  seismograms  for  (a)  model  2a,  (b)  model 
2b,  (c)  model  2c,  and  (d)  model  2d  (see  Table  1).  Solid  and  dotted  curves  stand 
for  the  phase  screen  and  finite  difference  solutions,  respectively. 


1162 


Y.-B.  Liu  and  R.-S.  Wu 


random  media,  (4)  exponential  random  media,  and  (5) 
flicker-noise  random  media,  which  were  called  “self- 
similar  random  media”  by  Frankel  and  Clayton  (1986). 
Figure  1  shows  typical  realizations  of  the  velocity  struc¬ 
ture  of  the  five  kinds  of  models.  There  are  considerably 
more  small-scale  structures  in  flicker-noise  media  than 
in  exponential  and  Gaussian  media.  For  the  definition 
and  generation  of  those  three  kinds  of  random  media  in 
the  2D  case,  see  Frankel  and  Clayton  (1986). 

For  absorbing  media,  we  assume  AzimFs  attenua¬ 
tion  law  (Azimi  et  al 1968;  Aki  and  Richards,  1980) 
is  satisfied  to  keep  the  propagation  causal.  The  corre¬ 
sponding  complex  wavenumber  and  complex  velocity 
have  the  following  dispersion  relation: 


where  a0  and  a,  are  constants  determining  the  absorbing 
properties  of  the  media,  and  v(°°)  is  the  limit  of  v(o>)  as 
w  oo  q  has  an  approximate  value  [2v(°°)a0]_1. 

In  order  to  reduce  the  boundary  effects,  we  choose 
the  model  space  long  enough  in  the  horizontal  direction 
and  make  observations  only  in  the  middle  section.  The 
model  sizes  are  51.2  by  51.2  km2  for  the  absorbing  cyl¬ 
inder  model  and  102.4  by  51.2  km2  for  the  other  four 
models.  Arrays  of  receivers  were  located  in  the  middle 
of  the  surface  along  the  x  direction.  A  schematic  dia¬ 
gram  of  the  model  space  is  shown  in  Figure  2.  The  spac¬ 
ing  of  the  grid  points  is  dx  =  dz  =  0.2  km. 

Table  1  lists  the  parameters  of  the  models  in  this 
study.  These  models  are  chosen  to  expose  how  the  ac¬ 
curacy  depends  on  the  correlation  function  (or  power 
spectrum),  magnitude  of  the  velocity  perturbation,  and 
the  ratio  of  wavelength  to  correlation  length  (or  scale 
length).  The  background  velocity  used  throughout  this 
study  is  v0  =  6.0  km /sec.  Correlation  lengths  of  5  km 
and  2  km  are  used.  The  values  k^a  are  given  in  Table 
1 ,  where  Jcq  denotes  the  wavenumber  of  dominant  fre¬ 
quency  (2  Hz). 


Distance(km) 


Distance(km) 


(a)  Gaussian  Random  Media  (e  =0.05,a=5km) 


(c)  Gaussian  Random  Media  (e  =0.15, a=5km) 


Distance  (km)  Distancc(km) 


(b)  Gaussian  Random  Media  (  E=0. 1  ,a=5km)  (d)  Gaussian  Random  Media  (  e=0.05,a=2km) 

Figure  6.  Comparison  of  synthetic  seismograms  for  (a)  model  3a,  (b)  model 
3b,  (c)  model  3c,  and  (d)  model  3d  (see  Table  1).  Solid  and  dotted  curves  stand 
for  the  phase  screen  and  finite  difference  solutions,  respectively. 
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. .  Source  Time  Function  and  Wavenumber  Filtering 

The  source  time  function  in  our  study  is  a  Ricker 
wavelet  defined  by 

g{t)  =  e-|/v*““<'-'0)2  cos  [n fndf  -  /0)],  (15) 

where is  the  maximum  frequency  and  t0  defines  the 
origin  of  time.  The  value  fmax  =  4  Hz  was  used  for  all 
simulations,  corresponding  to  a  central  frequency  of/0 
=  2  Hz.  The  wave  is  incident  from  the  bottom  of  the 
model.  The  same  time  step  A t  =  0.02  was  used  for  the 
phase  screen,  finite  difference,  and  eigenfunction  ex¬ 
pansion  methods. 

The  phase  screens  for  all  of  the  models  are  uni¬ 
formly  spaced  every  16  grid  points  (Az  =  3.2  km)  in  the 
z  direction.  This  spacing  is  used  in  order  to  satisfied  the 
geometrical  optics  condition  -within  screen  interval. 
Comparison  between  phase  screen  seismograms  with 
different  screen  intervals  has  been  done  and  showed  that 
there  was  no  noticeable  difference  between  seismograms 
of  up  to  16  grid  point  screen  intervals  (see  Figs.  3c  and 
3d).  In  order  to  reduce  the  artifacts  due  to  the  aliasing 


Distance  (km) 


(a)  Exponential  Random  Media  (e  =0.05^=5km) 


Distance(km) 


23  26  29  32  35  38  <1  44  47  50  53  56  39  62  «  «  71  74  77  80 


(b)  Exponential  Random  Media  (  e=0.1,a=5km) 


effect  at  very  high  wavenumbers,  a  wavenumber  filter 
is  applied  to  the  wave  field  in  wavenumber  domain  at 
each  step  (Nautiyal,  1988;  Wu  and  Huang,  1992).  Fig¬ 
ures  3a  and  3b  show  the  synthetic  seismograms  without 
wavenumber  filters  for  the  multi-uniform-cylinder  model 
and  the  flicker-noise  random  medium,  respectively.  In 
the  figures,  the  wiggles  before  the  first  arrivals  and  in 
the  later  codas  are  the  aliasing  artifacts.  After  applying 
the  wavenumber  filtering,  the  artifacts  are  effectively  re¬ 
moved  (Figs.  3c  and  3d).  Space-domain  cosine  tapers 
are  also  applied  to  the  wave  field  in  space  domain  near 
the  two  side  walls  to  reduce  the  side-wall  boundary  ef¬ 
fect. 

Comparison 

Since  the  finite  difference  and  eigenfunction  expan¬ 
sion  results  were  used  as  the  references  in  this  compar¬ 
ison,  it  is  important  to  assess  the  accuracy  of  those  meth¬ 
ods.  The  eigenfunction  expansion  is  an  analytic  solution 
and  can  give  exact  results  as  long  as  enough  terms  are 
taken.  To  test  the  degree  to  which  FD  converges  to  the 


Distance  (Von) 


(c)  Exponential  Random  Media  (E  =0.15,a=:51un) 


DistanceOcm) 


(d)  Exponential  Random  Media  (e  =0.05,a=2km) 


Figure  7.  Comparison  of  synthetic  seismograms  for  (a)  model  4a,  (b)  model 
4b,  (c)  model  4c,  and  (d)  model  4d  (see  Table  1).  Solid  and  dotted  curves  stand 
for  the  phase  screen  and  finite  difference  solutions,  respectively. 


1164 


Y.-B.  Liu  and  R.-S.  Wu 


exact  solution,  we  doubled  the  number  of  grid  points  for 
one  of  the  models  used  in  the  FD  calculations  and  found 
that  the  differences  in  the  solutions  were  negligible. 

Comparison  of  Phase  Screen  with  Finite 

Difference  for  Multi-Cylinder  Models 

Figure  4  gives  the  gray  scale  graphics  of  synthetic 
seismic  sections  from  finite  difference  simulation  and  the 
phase  screen  method  for  a  multi-uniform-cylinder  model 
(Figs.  4a  and  4b)  and  for  flicker-noise  random  media 
(Figs.  4c  and  4d),  respectively.  We  see  that  for  large- 
scale  heterogeneities,  such  as  those  in  the  multi-cylinder 
model,  the  agreement  between  phase  screen  and  FD  is 
excellent,  because  forward  scattering  is  dominant  in  this 
case.  On  the  other  hand,  the  agreement  deteriorates  for 
the  case  of  the  flicker-noise  medium,  especially  for  later 
arrivals,  since  the  flicker-noise  medium  is  rich  in  small- 
scale  heterogeneities  that  cause  strong  large-angle  scat¬ 
tering.  It  can  be  also  seen  that  because  of  the  treatment 
at  the  boundaries  (tapering,  absorbing  boundary  condi¬ 
tion),  the  amplitude  of  the  field  is  distorted  near  the 
boundaries.  Therefore,  in  the  following  we  will  only 


compare  the  synthetic  seismograms  away  from  the  ver¬ 
tical  boundaries;  i.e.,  in  the  regions  of  12  to  40  km  for 
the  absorbing  cylinder  model  and  23  to  80  km  for  the 
other  four  models.  In  this  way  we  can  also  avoid  the 
wrap-around  effects  due  to  the  equivalent  periodic 
boundary  condition  of  the  phase  screen  method  and  the 
boundary  reflections  of  the  FD  method  due  to  the  im¬ 
perfect  absorbing  boundary  conditions. 

Figure  5  shows  synthetic  seismograms  of  high-  and 
low-velocity  multi-uniform-cylinder  models  calculated 
by  phase  screen  and  finite  difference  for  receivers  lo¬ 
cated  in  the  middle  of  the  section.  Solid  and  dotted  curves 
stand  for  the  phase  screen  and  finite  difference  solutions, 
respectively.  It  can  be  seen  that  for  velocity  perturba¬ 
tions  with  e  =  +30%  (Fig.  5a),  e  —  —30%  (Fig.  5b), 
e  =  +50%  (Fig.  5c),  and  6  =  -50%  (Fig.  5d),  the  phase 
screen  results  are  in  excellent  agreement  with  the  finite 
difference  solutions,  demonstrating  that  the  phase  screen 
method  can  give  good  results  for  strong  discrete  scat¬ 
tered.  However,  as  can  be  noticed  from  Fig.  5d,  the 
strong  later  arrivals  for  the  case  of  low-velocity  cylinders 
have  large  errors  compared  with  the  FD  solution. 


Distance(km) 


23  26  25  32  35  3«  41  *4  <7  50  53  56  59  62  65  «  71  74  77  K) 


Distance  (km) 


(a)  Flicker-Noise  Random  Media  (e  =0.05,a=5km) 


(c)  Flicker-Noise  Random  Media  (£=0.15,a=5km) 
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- Comparison  of  Phase  Screen  with  Finite 

Difference  for  Random  Medium  Models 

Synthetic  seismograms  from  phase  screen  and  finite 
difference  methods  for  Gaussian  random  media,  expo¬ 
nential  random  media,  and  flicker-noise  random  media 
are  shown  in  Figures  6  through  8.  Solid  and  dotted  curves 
stand  for  the  phase  screen  and  finite  difference  results, 
respectively.  Figures  6a,  6d,  7a,  7d,  8a,  and  8d  are  the 
seismograms  of  the  three  kinds  of  random  media  with 
velocity  perturbations  of  e  =  5%  and  scale  lengths  a  = 
5  and  2  km,  respectively.  The  comparison  shows  ex¬ 
cellent  agreement  in  the  case  of  weak  random  media  (e 
^  5%). 

Figures  6b,  7b,  and  8b  are  the  seismograms  of  the 
three  kinds  of  random  media  with  velocity  perturbation 


Figure  9.  Root  mean  square  differences  of  the 
phase  screen  and  finite  difference  seismograms  for 
the  nine  models.  The  rms  differences  are  nor¬ 
malized  by  the  maximum  amplitude  at  the  cor¬ 
responding  time  for  each  curve. 


e  =  10%,  respectively.  It  is  obvious  that  the  agreement 
between  the  waveforms  at  early  times  is  very  good,  but 
it  degrades  at  later  times.  It  is  apparent  that  the  coda  at 
later  times  for  these  models  come  from  large-angle  mul¬ 
tiple  scattering. 

Figures  6c,  7c,  and  8c  are  the  seismograms  of  the 
three  kinds  of  random  media  with  velocity  perturbation 
e  =  15%,  respectively.  The  agreement  between  phase 
screen  and  finite  difference  solutions  degenerates  more 
than  that  in  Figures  6b,  7b,  and  8b.  The  coda  of  the 
phase  screen  calculations  is  noticeably  less  than  that  of 
the  finite  difference  results.  Clearly,  large-angle  scatter¬ 
ing  and  backscattering  play  significant  roles  in  produc¬ 
ing  the  larger  amplitude  of  codas  in  the  FD  results. 

To  quantify  the  comparison  between  the  phase  screen 
and  finite  difference  results.  The  rms  differences  be¬ 
tween  the  phase  screen  and  finite  difference  seismo¬ 
grams  were  calculated.  Figures  9a  (e  =  0.05,  a  =  5  km), 
9b  (e  =  0. 10,  a  =  5  km),  and  9c  (e  =  0. 15,  a  =  5  km) 
show  the  relative  rms  differences  of  seismograms  as 
functions  of  time  for  nine  models.  The  rms  averages  are 
taken  over  the  282  receivers  and  normalized  by  the  max¬ 
imum  amplitude  of  the  seismograms  at  the  correspond- 

Distancc  (km) 


(a)  Cylinder  Model  cv0  *  6-Ofcm/i.v^  *  4.8  knV*.Q  «  20.0) 


Distance  (km) 


(b)  Absorbing  Cylinder  {v0  «  6.0km/i.  =  4.8  km/s.o*  «  5.2x10  J.a,  =  8.0x10  Q  *  20.0) 

Figure  10.  Comparison  of  synthetic  seismo¬ 
grams  for  (a)  model  la  with  (dotted)  and  without 
(solid)  intrinsic  attenuation,  (b)  model  lb  w-ith 
(dotted)  and  without  (solid)  velocity  dispersion. 
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ing  time.  It  is  apparent  that  weak  random  media  (c  ^ 
10%)  are  treated  more  accurately  by  the  phase  screen 
method  than  for  the  strong  random  media  (e  =  15%), 
and  the  seismograms  for  earlier  arrivals  (from  8  to  10 
sec)  compare  with  the  FD  solutions  more  favorably  for 
Gaussian  and  exponential  than  for  flicker-noise  random 
media.  Large-angle  scattering  and  backscattering  are  more 
significant  in  the  latter  cases. 

Computation  Speed 

For  the  numerical  examples  calculated  in  this  arti¬ 
cle,  the  grid  size  is  512  by  256.  In  a  SUN  IV  station  the 
CPU  time  for  the  2D  FD  is  7108  sec,  while  for  the  phase 
screen  method  with  screen  interval  of  16  grid  points,  the 
CPU  time  is  478  sec,  about  15  times  faster  than  the  FD 
algorithm.  We  have  run  several  examples  of  1024  by 
512  grid  size  in  a  SUN  SPARC  II  station.  The  CPU  times 
for  FD  and  phase  screen  (with  screen  interval  of  32  points) 
are  20,917  sec  and  367  sec,  respectively.  The  time  sav¬ 
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ings  is  about  57  times.  For  large  3D  problems,  the  sav¬ 
ing  in  CPU  time  and  computer  storage  is  expected  to  be 
much  greater. 

Comparison  of  Phase  Screen  with  Exact  Solutions 

for  an  Absorbing  Cylinder 

Figure  10a  shows  the  synthetic  seismograms  from 
the  exact  solution  for  a  single  cylinder  model  with  and 
without  intrinsic  attenuation.  Solid  and  dotted  curves  stand 
for  the  nonabsorbing  and  absorbing  cylinder  scatterers, 
respectively.  Due  to  the  absorption  effect,  the  ampli¬ 
tudes  from  the  absorbing  cylinder  are  less  than  those  from 
the  nonabsorbing  cylinder. 

In  order  to  have  a  causal  attenuation  signal,  a  small 
amount  of  dispersion  is  introduced  according  to  the  Kra¬ 
mer  s-Kronig  relation  (Aki  and  Richards,  1980).  Figure 
10b  shows  the  synthetic  seismograms  from  exact  solu¬ 
tion  for  the  absorbing  cylinder  model  with  and  without 
phase  velocity  dispersion.  Solid  and  dotted  curves  stand 
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(b)  Absorbing  Cylinder  <v„  *  6.0  km v*  =  7.8  kmfr.Oo  =  3.21*10 ■  \a,  -  2.0*10  \Q  =  20.0>  (d)  Absorbing  Cylinder  (v0  *  6.0  km/,,  v*  *  9.0  km/s.a,  =  2.78*10 •  \a,  -  20*10  ’ \  Q  =  20  0, 

Figure  1 1 .  Comparison  of  PHS  (phase  screen)  and  eigenfunction  expansion 
results  for  (a)  model  lc  with  a0  =  5.95  x  10-3  and  a{  =  2.0  x  10  4;  (b)  model 
Id  with  a0  =  3.21  X  10-3  and  a,  =  2.0  x  I0~4;  (c)  model  le  with  aQ  =  8.33 
X  10"3  and  a,  =  2.0  X  10"4;  and  (d)  model  If  with  a0  =  2.78  X  10“3  and  a, 

=  2.0  X  10“4.  Solid  and  dotted  curves  stand  for  the  phase  screen  and  exact 
solutions,  respectively. 
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for  the  nondispersive  and  dispersive  cylinders,  respec¬ 
tively.  It  can  be  seen  that  some  waveform  change  occurs 
due  to  the  dispersion  effect. 

Figures  11a  and  1  lb  show  the  seismograms  from  the 
absorbing  cylinders  (with  causal  dispersion,  Q  =  20.0, 
Pi  “  Pi)  with  velocity  perturbations  of  e  =  -30%  ( aQ 
=  5.95  X  1(T3,  a,  =  2.0  X  10~4)  and  +30%  (a0  =  3.2 
X  10-3,  a|  =  2.0  X  10-4),  respectively.  Solid  and  dotted 
curves  stand  for  the  phase  screen  and  exact  solutions, 
respectively.  The  phase  screen  results  are  in  good  agree¬ 
ment  with  the  exact  solutions. 

Figures  11c  and  1  Id  show  the  seismograms  from  the 
absorbing  cylinders  (with  causal  dispersion,  Q  =  20.0, 
pi  =  p2 )  with  velocity  perturbations  e  =  —50%  (ot0  = 
8.33  X  1(T\  a  =  2.0  x  10^)  and  +50%  (a0  =  2.78 
X  10~3,  ax  =  2.0  X  10-4).  Solid  and  dotted  curves  stand 
for  the  phase  screen  and  exact  solutions,  respectively. 
The  agreement  between  the  phase  screen  and  exact  so¬ 
lutions  is  better  for  the  high-velocity  cylinder  than  for 
the  low-velocity  cylinder. 

Conclusions 

From  the  comparison  of  the  phase  screen  method 
with  the  exact  solution  for  a  cylinder  and  the  FD  solu¬ 
tions  for  various  heterogeneous  media  we  can  draw  the 
following  conclusions.  The  phase  screen  solution  can  give 
good  results  for  strong  discrete  heterogeneous  media  with 
and  without  intrinsic  attenuation  if  the  size  of  the  scat¬ 
tered  are  large  compared  with  the  wavelength.  For  con¬ 
tinuous  random  media,  the  comparison  of  synthetic  seis¬ 
mograms  by  phase  screen  and  finite  difference  methods 
was  made  for  three  kinds  of  models:  (1)  Gaussian  ran¬ 
dom  media,  (2)  exponential  random  media,  and  (3)  flicker- 
noise  random  media.  Results  show  good  agreement  for 
weak  random  media  (velocity  perturbations  ^=10%), 
similar  to  the  results  obtained  by  Fisk  et  al.  (1992).  The 
phase  screen  method  is  most  accurate  for  the  early  part 
of  the  seismograms,  where  the  arrivals  have  undergone 
only  small-angle  scattering.  At  later  times,  especially  for 
the  case  of  strong  scattering,  backscattering  becomes  more 
important  and  the  phase  screen  method  becomes  less  ac¬ 
curate  compared  with  the  finite  difference  method  be¬ 
cause  of  the  neglect  of  backscattered  waves  in  the  phase 
screen  formulation.  Therefore,  wherever  the  backscat¬ 
tered  waves  are  insignificant  or  can  be  separated  from 
the  major  wave  groups  of  interest,  the  phase  screen  method 
can  be  used  to  generate  synthetic  seismograms  with  much 
faster  computation  speed  than  FD  calculations.  The  other 
advantage  of  the  phase  screen  method  is  the  capability 
of  modeling  arbitrarily  absorbing  inhomogeneous  media. 
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Effects  of  Crustal  Structure  under  the  Barents  and  Kara  Seas  on 
Short-Period  Regional  Wave  Propagation  for  Novaya  Zemlya 
Explosions:  Empirical  Relations 

by  Tianrun  Zhang  and  Thome  Lay 


Abstract  Short-period  seismic  recordings  at  regional  and  upper  mantle  dis¬ 
tances  from  underground  explosions  at  Novaya  Zemlya  demonstrate  that  prop¬ 
agation  across  the  continental  shelf  under  the  Barents  and  Kara  Seas  appears  to 
modify  the  partitioning  of  energy  between  Lg  and  Sn  phases  relative  to  purely 
continental  paths  in  the  Eurasian  crust.  While  the  underwater  segments  of  the 
paths  are  relatively  short,  variations  in  bathymetric  characteristics  from  path  to 
path  influence  the  regional  wave  field,  with  systematic  behavior  that  can  be 
used  to  establish  empirical  amplitude  corrections  for  regional  phases.  We  ana¬ 
lyze  a  large  set  of  Eurasian  recordings  to  explore  the  relationship  between  re¬ 
gional  phase  energy  partitioning  and  bathymetric  characteristics.  Maximum  water 
depth  along  the  path  is  the  most  influential  factor  for  the  Novaya  Zemlya  data. 
It  has  strong  linear  correlations  with  the  logarithmic  rms  amplitude  of  Lg  and 
the  ratios  Sn/Lg  and  P/Lg.  The  maximum  water  depth  probably  reflects  the  ex¬ 
tent  of  necking  of  the  crustal  wave  guide  under  the  continental  margin,  which 
may  disrupt  Lg  modes  resulting  in  Lg  to  Sn  scattering,  but  there  is  surprising 
sensitivity  to  small  variations  in  bathymetry.  Empirical  relations  like  those  found 
here  may  be  useful  for  nuclear  yield  estimation  and  discrimination  for  regions 
such  as  the  Korean  Peninsula  and  Persian  Gulf,  where  many  seismic  phases 
traverse  water-covered  continental  shelf  with  poorly  known  crustal  structure. 


Introduction 

The  short-period  regional  seismic  phase  Lg ,  com¬ 
prised  of  Rayleigh  wave  overtones  or  postcritical  mul¬ 
tiply  reflected  S  waves,  has  long  been  known  as  a  “con¬ 
tinental  phase.”  In  fact,  its  first  applications  were  for 
discriminating  continental  and  oceanic  crust  based  on  the 
presence  or  absence  of  the  Lg  phase,  respectively  (Press 
and  Ewing,  1952).  Oliver  et  al.  (1955)  used  such  bi- 
modal  classifications  to  investigate  the  crustal  structure 
in  the  Arctic,  while  Savarensky  and  Valdner  (1960) 
studied  the  Black  Sea  region.  These  studies  concluded 
that  in  no  case  does  Lg  propagate  through  crust  overlain 
for  any  significant  distance  by  water  deeper  than  1000 
fathoms  (1.8  km).  Thin  oceanic  crust  is  now  well  known 
to  inhibit  Lg  propagation  (e.g.,  Kennett  and  Mykkeltveit, 
1984),  but  the  effects  of  transitional  crustal  structures  are 
not  as  clear. 

It  was  quickly  recognized  that  propagation  across  a 
marginal  sea  or  continental  shelf  does  not  completely 
quench  Lg ,  but  can  reduce  its  amplitude.  A  number  of 
studies  have  been  conducted  for  Lg  traversing  paths  un¬ 
der  shallow  water.  Wetmiller  (1974)  investigated  the 


crustal  structure  in  the  Baffin  Bay  area,  and  Gregersen 
(1984)  studied  the  crustal  structure  near  Denmark  and 
the  North  Sea.  Kennett  et  al.  (1985)  then  examined  het¬ 
erogeneity  of  the  North  Sea  basin  using  variability  of  Lg 
transmission.  Recently,  Baumgardt  (1990)  explored  Lg 
blockage  and  scattering  by  the  Barents  Sea  and  White 
Sea,  recognizing  the  importance  of  energy  partitioning 
between  Lg  and  Sn  in  the  regional  wave  field.  He  inferred 
that  the  presence  of  thick  sedimentary  basins  plays  a  crit¬ 
ical  role  in  Lg  blockage.  Nuttli  (1988)  analyzed  Lg  phases 
on  similar  paths  across  the  Barents  Sea,  with  only  mod¬ 
erately  strong  attenuation  values  being  found  for  the 
blocked  arrivals,  and  their  relative  yield-amplitude  scal¬ 
ing  appearing  to  be  unimpaired.  Theoretical  studies  of 
propagation  across  continental  margins  have  been  con¬ 
ducted  by  Kennett  (1986),  Maupin  (1989),  Regan  and 
Harkrider  (1989),  and  Cao  and  Muirhead  (1993).  The 
latter  study  used  P-SV  finite  difference  calculations  to 
demonstrate  that  necking  of  the  crust  and  the  presence 
of  the  thin  water  layer  can  effectively  block  Lg.  Many 
additional  studies  have  examined  variations  in  Lg  and  Sn 
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transmission  efficiency  in  diverse  continental  environ¬ 
ments  (e.g.,  Ruzaikin  et  al. ,  1977;  Kadinsky-Cade  et 
al.,  1981;  Ni  and  Barazangi,  1983;  Lynnes  and  Baum- 
stark,  1991). 

These  important  studies  have  provided  insight  into 
the  complex  structural  controls  on  Lg  in  a  heterogeneous 
transitional  crust,  but  we  still  lack  a  practical  means  by 
which  to  account  for  the  strong  amplitude  and  spectral 
effects  of  underwater  segments  for  regional  arrivals.  This 
is  particularly  important  for  efforts  to  monitor  under¬ 
ground  explosions  of  low  yield  and  to  discriminate  them 
from  earthquakes  using  regional  phases.  The  strong  vari¬ 
ability  of  regional  discriminants  (e.g.,  Pomeroy  et  al ., 
1982;  Lynnes  and  Baumstark,  1991)  mandates  correc¬ 
tion  for  path  properties  affecting  regional  phase  ampli¬ 
tudes  and  spectral  content. 

The  water  in  most  parts  of  the  North  Sea,  Barents 
Sea,  and  Kara  Sea  is  shallower  than  200  m.  The  thinning 
of  the  crustal  wave  guide  under  the  shelf  is  small  (except 
in  the  central  rift  of  the  North  Sea)  relative  to  the  typical 
30-  to  40-km  thickness  of  the  adjacent  continental  crust. 
However,  Lg  waves  do  have  highly  variable  transmission 
efficiency  through  these  water-covered  regions.  This  im¬ 
plies  that  there  must  be  significant  changes  in  the  crustal 
structure  or  in  the  shallow  sediments  underlying  the  mar¬ 
ginal  seas  that  affect  the  regional  wave  field. 

This  study  explores  the  influence  of  the  bathymetric 
structure  of  the  Barents  Sea  and  Kara  Sea  on  regional 
phase  propagation.  The  approach  taken  is  an  empirical 
one,  given  that  current  capabilities  for  realistically  sim¬ 
ulating  regional  waveforms  are  quite  limited,  as  is  our 
detailed  knowledge  of  the  regional  crustal  structure.  One 
of  the  major  obstacles  confronting  nuclear  nonprolifer¬ 
ation  monitoring  is  the  variability  of  regional  wave  en¬ 
ergy  flux  associated  with  diverse  paths  in  the  crust.  If  it 
is  possible  to  develop  reliable  empirical  approaches  that 
use  surface  observations  such  as  path  topography  or  ba¬ 
thymetry  to  account  for  some  of  the  wave  field  vari¬ 
ability,  more  effective  procedures  for  global  monitoring 
may  be  achieved.  For  example,  Zhang  and  Lay  (1994) 
have  recently  discovered  surprisingly  strong  empirical 
relationships  between  regional  phase-amplitude  ratios  and 
along-path  topography  statistics  for  Eurasian  continental 
paths.  In  this  spirit,  our  analysis  involves  seeking  simple 
path  parameters  that  provide  a  basis  for  reducing  wave 
field  variations  associated  with  propagation  through  the 
crustal  structure  underlying  a  shallow  sea.  Quantitative 
modeling  of  any  such  relationships  will  be  pursued  in 
the  future. 

Data 

The  seismic  data  that  we  use  are  recordings  of  20 
underground  nuclear  explosions  located  in  the  Novaya 
Zemlya  test  site  near  Matochkin  Shar  (the  narrow  strait 
that  separates  the  northern  and  southern  islands  of  No¬ 


vaya  Zemlya).  The  events  occurred  from  1964  to  1988. 
Short-period  regional  and  upper  mantle  distance  record¬ 
ings  from  Soviet-run  stations  were  collected  and  digi¬ 
tized  for  these  events  as  part  of  a  data  exchange  asso¬ 
ciated  with  the  Joint  Verification  Experiment.  The 
locations  of  the  test  site  and  stations  providing  sufficient 
numbers  of  recordings  and  instrument  response  infor¬ 
mation  for  our  analysis  are  shown  in  Figure  1.  Topog¬ 
raphy  contours  for  Eurasia  are  included,  derived  from 
digital  elevation  data  supplied  at  5-min  intervals  from  a 
world  topography  data  base  (ET0P05,  compiled  by  the 
National  Geophysical  Data  Center,  Boulder,  Colorado). 
The  stations  range  in  distance  from  1900  to  3600  km, 
and  provide  the  most  extensive  data  set  for  Novaya  Zem¬ 
lya  explosions  available  from  the  Eurasian  mainland.  Most 
previous  studies  of  short-period  regional  phases  from  ex¬ 
plosions  at  this  test  site  have  used  array  data  from  Nor¬ 
way,  Germany,  and  Finland  (e.g.,  Baumgardt,  1990; 
Ringdal  and  Fyen,  1991)  or  from  isolated  stations  in 
western  Europe  and  Scandinavia  (Nuttli,  1988),  which 
sample  a  very  restricted  azimuthal  range  from  the  test 
site.  The  exchange  data  provide  an  extensive  data  set 
with  broad  azimuthal  coverage,  including  paths  with  di¬ 
verse  underwater  segments. 

Figure  2  shows  detailed  bathymetry  of  the  region 
around  Novaya  Zemlya,  with  the  Barents  Sea  to  the  west 
and  the  Kara  Sea  to  the  east.  The  ray  paths  to  the  ex¬ 
change  stations  sample  different  paths  through  these 
marginal  seas.  This  provides  a  basis  for  exploring  the 
relative  waveform  effects  of  different  underwater  seg¬ 
ments.  The  previous  work  on  Novaya  Zemlya  regional 
data  mentioned  above  has  primarily  used  paths  travers¬ 
ing  the  deeper  water  regions  of  the  Barents  Sea  (toward 
the  west  in  Fig.  2). 

Only  vertical-component  seismograms  have  been 
digitized  for  the  exchange  stations,  but  the  general  signal 
quality  is  quite  good  (Israelsson,  1992).  Altogether,  108 
waveforms  from  the  20  explosions  recorded  by  eight  sta¬ 
tions  are  used  in  the  analysis.  The  instrumentation  is 
comparable  between  stations,  but  there  are  slight  vari¬ 
ations  in  frequency  response,  as  well  as  variations  with 
time  at  each  station.  To  equalize  the  instrument  char¬ 
acteristics,  we  deconvolved  the  individual  instrument  re¬ 
sponses  from  each  waveform,  and  then  convolved  the 
ground  motion  with  the  1988  response  of  the  CKM-3 
seismometer  at  station  OBN.  For  some  of  the  data,  por¬ 
tions  of  the  records  are  missing;  we  use  only  the  phases 
which  are  reliably  digitized.  Many  details  about  the  spe¬ 
cific  explosions,  the  exchange  data  set,  and  the  signal 
quality  are  given  by  Israelsson  (1992).  We  combine  data 
from  different  events  to  characterize  average  path  prop¬ 
erties,  and  so  do  not  dwell  on  the  individual  event  in¬ 
formation. 

Our  primary  interest  is  to  explore  whether  energy 
partitioning  in  the  regional  wave  field  is  influenced  by 
the  variation  in  underwater  path  segments.  To  achieve 
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this,  we  require  robust  measurements  of  the  short-period 
wave  field.  Previous  studies  have  established  the  stabil¬ 
ity  of  gross  averaging  measures  such  as  rms  calculations 
for  regional  phases  (e.g.,  Ringdal  et  al.y  1992;  Hansen 
et  al. ,  1990;  Israelsson,  1992;  Zhang  and  Lay,  1994), 
so  we  use  comparable  measurements  in  this  study.  The 
rms  values  for  different  phases  were  calculated  using  the 
instrument-equalized  traces,  with  bandpass  filtering  be¬ 
tween  0.6  and  3.0  Hz.  These  frequency  limits  conform 
to  routine  regional  phase  passbands  used  by  Hansen  et 
al.  (1990)  and  are  compatible  with  the  limitations  of  the 
hand-digitized  data. 

The  rms  amplitudes  of  each  phase  are  calculated  in 
corresponding  group  velocity  or  time  windows.  The  Lg 
is  assigned  the  window  3.i  to  3.7  km/sec,  following 
Israelsson  (1992),  while  the  Sn  window  is  4.3  to  4.8  km/ 
sec,  as  preferred  by  Kennett  (1989)  (see  Fig.  4  for  an 
example).  We  found  that  use  of  a  narrower  Lg  window 
of  3.1  to  3.5  km/sec  had  only  minor  effects  on  the  mea¬ 
surements  in  this  study.  We  use  time  windows  for  the  P 
waves,  which  at  the  upper  mantle  distances  spanned  by 


our  data  have  impulsive  first  arrivals  followed  by  com¬ 
plex  coda.  We  use  both  0-  to  5-  and  0-  to  50-sec  win¬ 
dows  from  the  onset  of  the  P  arrivals.  The  latter  part  of 
the  0-  to  50-sec  window  contains  an  increasing  com¬ 
ponent  that  is  multiply  reflected  within  the  crustal  wave 
guide,  but  the  primary  energy  dives  into  the  upper  man¬ 
tle.  A  noise  correction  was  applied  to  each  measurement 
(Zhang  and  Lay,  1994),  obtained  from  the  rms  ampli¬ 
tude  of  the  available  signal  prior  to  5  sec  preceding  the 
manually  picked  P  onset  (usually  115  sec  of  data).  A 
comparably  measured  data  set  for  underground  explo¬ 
sions  at  the  Semipalatinsk  test  site,  comprised  of  325 
waveforms  for  83  explosions  from  the  JVE  exchange 
(Zhang  and  Lay,  1994)  is  used  for  comparison  with  the 
Novaya  Zemlya  observations.  The  paths  associated  with 
those  data  are  indicated  in  Figure  1 . 

Characterization  of  Paths 

The  islands  of  Novaya  Zemlya  are  an  extension  of 
the  Ural  mountains,  commonly  identified  with  the 


Figure  1 .  Map  showing  the  locations  of  the  Novaya  Zemlya  and  Semipala¬ 
tinsk  test  sites,  marked  with  triangles,  and  the  seismological  stations  used,  marked 
with  circles  and  indicated  by  their  codes.  The  lines  between  the  sources  and 
receivers  correspond  to  great-circle  paths.  Topography  is  indicated  with  contour 
intervals  every  1000  m  from  -1000  to  3000  m  above  sea  level,  but  500  m  con¬ 
tour  is  added  to  show  Ural  and  other  low  mountains.  Plotting  software  from 
Wessel  and  Smith  (1991)  is  used  in  this  and  the  next  two  figures. 
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Figure  2.  Map  showing  bathymetry  around  Novaya  Zemlya,  and  great-circle 
paths  from  the  test  site  to  each  Eurasian  station  used.  The  contour  interval  as¬ 
signed  to  each  gray-shade  is  50  m.  The  Barents  Sea  is  to  the  west  of  the  island, 
and  the  Kara  Sea  to  the  east.  The  paths  traverse  the  continental  shelf  with  depths 
ranging  from  50  m  to  a  maximum  of  350  m  (east  of  the  island). 


Table  1 

Bathymetric  Statistics  of  the  Paths 


Path 

Mean  Depth  (m) 

Length  (km)* 

Max.  Depth  (m) 

Area  (km')’ 

BOD 

111 

420 

341 

46.9 

NRI 

108 

420 

338 

45.3 

TLY 

126 

420 

278 

52.8 

NVS 

116 

590 

337 

68.4 

FRU 

220 

280 

340 

61.6 

ARU 

41 

230 

116 

9.5 

OBN 

6 

450 

129 

27.1 

UZH 

87 

710 

171 

61.7 

*Total  underwater  length. 

tCross-sectional  area  of  the  underwater  segment. 


boundary  between  Europe  and  Asia.  The  eight  paths  from 
Novaya  Zemlya  to  Eurasian  stations  sampled  by  our  data 
are  readily  divided  into  two  groups.  Three  paths,  to  UZH, 
OBN,  and  ARU,  cross  the  southeast  Barents  Sea  (Fig. 
2),  where  the  maximum  water  depth  is  100  to  150  m. 
We  call  this  the  Barents  group.  There  is  a  relatively  deep¬ 


water  zone  off  the  east  coast  of  Novaya  Zemlya,  in  the 
southern  Kara  Sea.  Paths  to  FRU,  NVS,  TLY,  NRI,  and 
BOD  traverse  this  deeper  water,  where  the  maximum  depth 
is  about  300  m  (Fig.  2).  This  is  the  Kara  group.  Each 
path  has  slightly  different  underwater  segment  charac¬ 
teristics,  which  we  will  consider.  These  include  mean 
depth,  total  underwater  length,  maximum  depth,  and  the 
cross-sectional  area  of  the  underwater  segment,  the  val¬ 
ues  for  which  are  given  in  Table  1. 

The  entire  propagation  path  influences  the  regional 
wave  field,  and  in  this  case  all  of  the  paths  have  signif¬ 
icant  segments  across  the  Russian  platform.  While  de¬ 
tails  of  the  crustal  structure  on  each  path  are  at  best  very 
sketchy,  a  possible  guide  to  relative  structural  charac¬ 
teristics  is  provided  by  surface  topography.  The  topog¬ 
raphy  profiles,  from  the  ETOP05  data  base  with  10-km 
lateral  sampling,  are  shown  for  each  path  in  Figure  3. 
The  island  source  region  is  on  the  left  of  each  profile 
and  the  station  is  on  the  right.  There  is  a  large  vertical 
exaggeration,  but  in  most  cases  the  underwater  segments 
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are  minor  features  relative  to  the  continental  relief  along 
the  paths. 

The  path  to  BOD  has  a  mean  water  depth  of  1 1 1  m, 
and  a  maximum  water  depth  of  341  m,  similar  to  the 
other  four  paths  in  the  Kara  group.  The  pathlength  across 
the  sea  is  420  km.  We  do  not  include  the  pathlength 
under  shallow  water  bodies  along  this  path  like  Obskaya 
bay  and  Yeniseyskiy  gulf.  The  western  half  of  the  BOD 
on-land  path  is  flat,  on  the  Siberian  lowland.  When  the 
path  intersects  the  Vilyuyskoye  Plateau,  the  altitude  rises 
and  the  surface  becomes  quite  rough.  On  the  basis  of  the 
work  by  Zhang  and  Lay  (1994),  we  anticipate  that  this 
topography  will  reflect  changes  in  the  wave  guide  that 
affect  the  propagation  of  regional  phases,  although  pos¬ 
sibly  not  as  dramatically  as  the  underwater  segment.  The 
path  to  NRI  almost  coincides  with  the  west  half  of  the 
BOD  path,  ending  just  before  the  rough  topography  be¬ 
gins.  The  path  to  TLY  is  similar  to  that  to  BOD,  with 
the  southern  segment  encountering  rough  topography  in 
the  Sayan  mountains.  The  path  to  NVS  has  a  long  un~ 
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Figure  3.  Profiles  of  surface  topography  from 
the  ET0P05  data  base  along  paths  between  the 
Novaya  Zemlya  test  site  and  the  eight  Eurasian 
stations  shown  in  Figure  1.  The  units  of  horizontal 
and  vertical  axes  are  both  kilometers,  but  the  ver¬ 
tical  scale  is  exaggerated  by  a  factor  of  110.  The 
dotted  lines  are  sea  level. 


derwater  segment,  as  it  traverses  Baydaratslaya  bay,  the 
southernmost  part  of  the  Kara  Sea.  The  on-land  segment 
is  quite  flat.  The  path  to  FRU  has  the  greatest  mean  water 
depth  (220  m)  and  a  maximum  depth  of  340  m  (Table 
1).  Its  on-land  segment  crosses  the  Ural  mountains  and 
then  traverses  the  huge  Siberian  plane,  before  encoun¬ 
tering  the  high-altitude  Pamir  mountains.  Unfortunately, 
there  is  only  one  record  available  from  FRU. 

The  paths  of  the  Barents  group  have  relatively  shal¬ 
low  underwater  segments.  The  path  to  ARU  runs  along 
the  Ural  mountains  for  a  few  hundred  kilometers,  which 
is  likely  to  affect  the  regional  phase  propagation.  The 
path  to  OBN  has  a  very  shallow  mean  water  depth  and 
a  flat  on-land  path.  The  path  to  UZH  has  the  longest  un¬ 
derwater  segment  (710  km),  but  the  mean  depth  is  only 
87  m.  This  path  crosses  the  White  Sea  and  Lake  Ladoga, 
but  we  have  only  summed  the  pathlength  under  the  Bar¬ 
ents  Sea. 


Examination  of  Path  Effects 

Oliver  et  al.  (1955)  observed  clear  Lg  arrivals  at  a 
station  in  Copenhagen  from  an  earthquake  located  in  the 
Arctic  Ocean,  with  a  path  that  crosses  the  Kara  Sea,  No¬ 
vaya  Zemlya  Island,  and  the  Barents  Sea.  Their  conclu¬ 
sion,  and  all  other  lines  of  evidence,  indicate  that  the 
crust  underlying  these  two  seas  is  essentially  continental. 
The  water  depth  between  Novaya  Zemlya  and  the  Eur¬ 
asian  mainland  is  never  more  than  350  m,  and  Lg  energy 
is  not  totally  blocked.  Yet,  we  find  that  traversing  the 
underwater  segment  appears  to  substantially  modify  the 
energy  partitioning  in  the  short-period  regional  wave  field, 
as  has  been  previously  demonstrated  for  paths  to  Europe 
(e.g.,  Baumgardt,  1990).  The  clearest  evidence  for  this 
is  provided  by  comparison  of  Novaya  Zemlya  recordings 
with  Eurasian  recordings  of  events  at  other  Soviet  test 
sites. 

In  order  to  evaluate  the  effects  of  the  underwater 
paths,  we  compare  the  seismograms  from  Novaya  Zem¬ 
lya  events  with  recordings  of  Semipalatinsk  explosions 
from  the  same  Eurasian  stations.  The  latter  involve  al¬ 
most  entirely  on-land  paths  (Fig.  1).  Figure  4  shows  rep¬ 
resentative  waveforms  at  similar  distances  for  recordings 
of  a  Novaya  Zemlya  event  at  station  OBN  (top)  and  a 
Semipalatinsk  event  at  the  same  station  (bottom).  To  em¬ 
phasize  the  relative  energy  distribution  in  the  short-pe¬ 
riod  signals,  the  times  scales  are  modified  to  have  com¬ 
mon  group  velocity  scales  from  9.0  to  2.8  km/sec.  The 
general  structural  characteristics  of  each  path,  as  mani¬ 
fested  in  surface  topography  and  geologic  provinces  tra¬ 
versed,  are  quite  similar,  apart  from  the  short  underwater 
segment  for  the  Novaya  Zemlya  data  (Fig.  1).  The 
waveform  from  the  Semipalatinsk  event  has  a  clear  Lg 
wave  packet  in  the  group  velocity  window  of  3.7  to  3.1 
km /sec,  with  a  maximum  near  3.5  km/sec.  The  Sn  win¬ 
dow  has  very  weak  energy,  indistinguishable  in  this 
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passband  from  the  P  coda.  The  Lg  window  for  the  No- 
vaya  Zemlya  event  has  lower  signal  to  noise  ratio  but 
there  is  clearly  energy  in  the  interval.  There  is  relatively 
high-amplitude  energy  in  the  latter  part  of  the  Sn  win¬ 
dow.  The  corresponding  path  traverses  the  shallowest 
water  segment  amongst  the  various  Novaya  Zemlya  paths, 
but  even  in  this  case  the  overall  difference  in  the  short- 
period  energy  distribution  is  apparent. 

The  differences  seen  in  Figure  4  are  common  to  the 
large  data  sets  for  each  test  site.  To  illustrate  this  we  use 
a  simple  measure,  picking  the  group  velocity  at  which 
each  waveform  amplitude  peaks  in  the  range  5.0  to  2.8 
km/sec.  Comparisons  of  such  measurements  for  a  large 
suite  of  explosions  at  the  two  sites  are  shown  in  Figure 
5.  Crosses  indicate  the  group  velocity  of  the  peak  rms 
amplitude  in  a  moving  8-sec  window  for  Semipalatinsk 
recordings,  while  circles  are  for  Novaya  Zemlya  data. 
Almost  all  Semipalatinsk  measurements  are  within  the 
3.7  to  3.1  km/sec  range,  corresponding  to  the  Lg  win¬ 
dow.  However,  most  Novaya  Zemlya  measurements  are 
in  the  4.6  to  3.8  km/sec  range,  which  is  the  latter  part 
of  the  Sn  (4.8  to  4.3  km/sec)  window,  and  in  what  Zhang 
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Figure  4.  Representative  seismograms  for  re¬ 
cordings  on  the  Semipalatinsk-OBN  path  {below) 
and  the  Novaya  Zemlya-OBN  path  {above).  The 
waveforms  are  bandpass-filtered  between  0.6  and 
3.0  Hz.  The  time  scales  are  slightly  different  so 
that  the  two  signals,  which  are  at  different  dis¬ 
tances  from  the  sources,  can  be  put  on  the  same 
group  velocity  scale  (indicated  at  the  top  of  each 
plot).  The  vertical  lines  mark  the  windows  used 
to  calculate  rms  amplitudes.  The  S  to  N  interval 
corresponds  to  the  Sn  window  of  4.8  to  4.3  km/ 
sec,  and  L  to  G  indicates  the  Lg  window  of  3.7 
to  3.1  km/sec.  Notice  the  relatively  strong  S„  for 
the  Novaya  Zemlya  event  recording. 


and  Lay  (1994)  define  as  the  Sn  coda  window  (4.3  to  3.7 
km/sec).  While  the  separation  is  on  average  very  clear, 
despite  the  variety  of  paths  involved,  there  are  a  few 
interesting  exceptions.  For  the  Semipalatinsk  data,  the 
two  records  on  the  path  to  APA  have  peak  amplitudes  at 
group  velocities  higher  than  4.0  km/sec.  This  path  is 
very  long  and  unusual  in  that  it  traverses  the  White  Sea 
just  before  reaching  the  station  (Fig.  1).  Although  the 
maximum  depth  of  the  White  Sea  is  only  48  m,  the  fact 
that  Sn  is  larger  than  Lgy  as  is  observed  for  Novaya  Zem¬ 
lya  records,  is  noteworthy.  For  the  Novaya  Zemlya  data, 
some  of  the  records  at  stations  ARU  and  FRU  have  peak 
amplitudes  at  low  group  velocities.  These  paths  are  close 
to  the  underwater  extension  of  the  Ural  mountains,  where 
the  crust  is  likely  to  be  relatively  thick. 

The  distance  ranges  spanned  by  the  Semipalatinsk 
and  Novaya  Zemlya  recordings  are  comparable,  but  there 
are  some  differences  in  the  amplitude  variations  of  the 
regional  phases  with  distance.  For  Novaya  Zemlya,  the 
rms  amplitude  in  the  Sn  window  of  4.8  to  4.3  km/sec  is 
considerably  enhanced  relative  to  Semipalatinsk  data.  In 
Figure  6,  squares  indicate  the  mean  rms  Sn  amplitudes 
from  Semipalatinsk,  while  circles  are  for  Novaya  Zem¬ 
lya.  The  values  are  computed  after  equalizing  the  event 
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Figure  5.  A  summary  plot  indicating  the  group 
velocity  at  which  a  running  average  amplitude  in 
an  8-sec  window  has  a  maximum  value  for  each 
seismogram.  The  measurements  are  plotted  versus 
distance  from  each  test  site,  with  circles  for  No¬ 
vaya  Zemlya  data  (stations  are  labeled  at  the  top), 
and  crosses  for  Semipalatinsk  data  (stations  are 
labeled  at  the  bottom).  The  group  velocity  range 
considered  is  from  5.0  to  3.0  km/sec. 
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sizes  by  scaling  the  measurements  to  a  common  teleseis- 
mic  mb  of  5.0.  The  Sn  is  generally  stronger  for  the  paths 
that  cross  underwater  segments.  Strong  Sn  was  observed 
by  Kadinsky-Cade  et  aL  (1981,  Figs.  12,  13  of  that  ar¬ 
ticle),  for  paths  across  the  Caspian  Sea  and  Black  Sea. 
They  found  that  the  attenuation  factors  for  Sn  and  Lg  are 
not  consistent  in  that  region,  but  they  did  not  conclude 
that  the  Sn  is  enhanced.  Baumgardt  (1990)  suggests  that 
the  late  Sn  and  Sn  coda  windows  may  be  enhanced  for 
paths  traversing  underwater  segments  due  to  Lg  to  Sn 
scattering,  so  it  is  possible  that  the  Novaya  Zemlya  rms 
values  are  enhanced  in  an  absolute  sense.  It  is  also  clear 
from  Figure  6  that  the  Sn  energy  for  Novaya  Zemlya  ex¬ 
hibits  a  different  distance  decay.  Some  of  the  Sn  obser¬ 
vations  from  Semipalatinsk  .tend  to  traverse  regions  with 
rough  topography  associated  with  the  Sayan  mountains 
and  lake  Baikal  (Zhang  and  Lay,  1994),  thus,  their  am¬ 
plitudes  and  distance  behavior  may  be  anomalous,  rel¬ 
ative  to  the  comparatively  flat  paths  for  the  Novaya 
Zemlya  data. 

Similar  comparisons  of  mean  rms  Lg  amplitudes  on 
paths  for  Novaya  Zemlya  (solid  circles)  and  Semipala- 
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Figure  6.  A  plot  of  the  mean  Sn  amplitudes  and 
their  standard  deviations  on  each  path  versus  path 
distance,  where  all  events  have  been  equalized  to 
a  common  mb  =  5.0.  Squares  are  for  Semipala¬ 
tinsk  data,  and  circles  are  for  Novaya  Zemlya.  The 
observed  Sn  amplitudes  from  Novaya  Zemlya  tend 
to  be  higher  and  show  a  stronger  distance  trend. 
Regressions  for  the  trend  with  distance  are  shown. 
For  each  regression,  CC  gives  the  correlation 
coefficient.  SIG  the  standard  deviation,  and  SLO 
the  regression  slope.  The  regression  curves  are  in¬ 
dicated  by  solid  lines.  Dashed  lines  mark  the  one- 
standard-deviation  limits,  ±cr. 


tinsk  (open  squares)  events  are  shown  in  Figure  7.  The 
events  are  again  equalized  to  a  common  mb  of  5.0.  The 
absolute  amplitude  levels  are  not  systematically  different 
overall,  in  contrast  to  the  Sn  measurements;  however,  the 
Novaya  Zemlya  data  exhibit  somewhat  greater  variabil¬ 
ity.  We  calculated  the  standard  deviation  of  the  distance 
trends  for  the  Lg  amplitudes  for  Novaya  Zemlya  and 
Semipalatinsk  data  separately,  finding  that  the  standard 
deviation  for  Novaya  Zemlya  data  is  about  twice  that  for 
Semipalatinsk  (0.29  and  0.14,  Fig.  7).  The  long  and  short 
dashed  lines  indicate  the  standard  deviations  for  the  two 
data  sets,  respectively.  Four  of  the  Novaya  Zemlya  data 
points,  for  stations  NRI,  NVS,  BOD,  and  TLY,  are  near 
the  lower  dashed  line.  These  all  belong  to  the  “Kara 
group” — the  deep-water  group,  and  it  would  appear  that 
on  average  the  Lg  amplitudes  are  reduced  for  these  paths 
relative  to  the  Semipalatinsk  data.  Of  the  four  Novaya 
Zemlya  points  near  the  upper  line,  three  (ARU,  OBN, 
and  UZH)  are  from  the  “Barents  group” — the  shallow- 
water  group.  The  only  exception  to  this  grouping  is  FRU, 
for  which  there  is  only  a  single  observation.  Uncertainty 
in  instrument  gains  may  contribute  to  some  of  the  scatter 
in  Figures  6  and  7,  since  some  stations  tend  to  plot  con¬ 
sistently  high  or  low  in  each  population. 

While  there  are  clearly  many  possible  factors  influ¬ 
encing  these  regional  phases,  we  are  very  limited  in  our 
ability  to  quantitatively  model  almost  all  of  them,  in¬ 
cluding  likely  factors  such  as  differences  in  excitation 
associated  with  the  near-source  properties  at  the  two  test 
sites.  We  will  emphasize  the  relative,  apparently  path- 
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Figure  7.  Similar  to  Figure  6,  but  for  Lg  am¬ 
plitudes.  The  observed  Lg  amplitudes  from  No¬ 
vaya  Zemlya  are  scattered  more  than  those  from 
Semipalatinsk. 
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dependent  variations,  as  these  can  be  treated  somewhat 
independently  of  the  absolute  excitation  issues.  In  pre¬ 
vious  studies  of  the  Semipalatinsk  data  set,  the  Lg  phase 
has  been  demonstrated  to  be  relatively  stable  (Israelsson, 
1992;  Zhang  and  Lay,  1994),  but  the  ratio  of  Sn/Lg  was 
found  to  have  a  strong  correlation  with  the  mean  altitude 
and  roughness  of  the  path  to  each  station  (Zhang  and 
Lay,  1994).  Figure  8  shows  the  correlation  of  the  rms 
amplitude  ratio  Sn/Lg  with  average  path  roughness  (stan¬ 
dard  deviation)  estimated  from  the  ET0P05  data  base, 
for  the  Semipalatinsk  and  Novaya  Zemlya  data  sepa¬ 
rately.  Note  the  very  strong  correlation  for  Semipala¬ 
tinsk,  and  the  relatively  large  range  of  roughness  sam¬ 
pled  by  the  different  paths  (details  are  given  in  Zhang 
and  Lay,  1994).  One  interpretation  is  that  the  surface 
roughness  provides  a  surrogate  for  irregularities  of  the 
crustal  wave  guide,  which  influence  the  energy  parti¬ 
tioning  between  Sn  and  Lg  along  each  path.  In  contrast, 
the  Novaya  Zemlya  data  do  not  show  a  corresponding 
relationship,  despite  the  fact  that  many  of  the  paths  tra¬ 
verse  similar  Eurasian  structure.  In  part,  this  may  reflect 
the  fact  that  the  overall  variation  in  path  roughness  is 
reduced  for  the  Novaya  Zemlya  data  set,  but  we  suspect 
that  a  more  important  factor  is  the  variability  associated 
with  the  differences  in  underwater  segments,  which  ap¬ 
pears  to  overwhelm  the  on-land  effects.  We  will  test 
whether  we  can  isolate  any  such  effect. 

We  use  the  four  underwater  path  parameters,  mean 
depth,  underwater  length,  maximum  depth,  and  under¬ 
water  cross-sectional  area,  tabulated  in  Table  1 ,  to  char¬ 
acterize  the  variations  in  shelf  segment.  Our  strategy  is 
to  correlate  the  amplitudes,  and  amplitude  ratios,  of  each 
phase  with  the  four  factors.  Figure  9  shows  the  result  for 
average  rms  Lg  measurements  at  each  station,  which  have 
been  normalized  to  a  common  mb  =  5.0  and  then  cor¬ 
rected  for  distance  according  to  power  law  r~s/6  (Nuttli, 
1973).  The  latter  correction  is  only  for  geometric  spread¬ 
ing,  and  ignores  effects  of  attenuation,  which  will  be 
considered  later.  Mean  depth  (Fig.  9a)  and  underwater 
area  (Fig.  9d)  have  moderate  correlation  with  Lg  ampli¬ 


tude  variations.  The  length  of  water  segment  (Fig.  9b) 
shows  no  pattern.  Maximum  water  depth  (Fig.  9c)  shows 
the  strongest  pattern  among  the  four  factors,  with  the 
Barents  group  (ARU,  OBN,  and  UZH)  clustering  on  the 
left,  while  the  Kara  group  (TLY,  NVS,  NRI,  FRU,  and 
BOD)  clusters  on  the  right.  The  Kara  group  involves  ob¬ 
servations  spanning  a  large  distance  range.  A  factor  of 
3  difference  in  amplitudes  is  involved,  corresponding  to 
the  large  scatter  in  Figure  7.  Mean  depth,  maximum  depth, 
and  cross-sectional  area  are  all  correlated,  but  it  appears 
that  maximum  depth  is  particularly  important,  more  so 
than  the  length  of  the  water  segment.  Mean  depth  would 
have  a  much  stronger  correlation  if  the  single  observa¬ 
tion  at  FRU  were  omitted. 

Very  similar  results  are  found  for  rms  Sn  amplitudes 
(Fig.  10).  The  same  geometric  spreading  correction  has 
been  applied  to  the  Sn  data.  The  Sn  amplitude  also  de¬ 
creases  with  increasing  maximum  water  depth,  but  the 
amplitude  variations  are  somewhat  less  than  for  Lg. 

If  wave  guide  properties  associated  with  water  depth 
are  in  fact  important  for  the  regional  wave  field,  we  would 
expect  there  to  be  no  correlation  with  rms  amplitude 
variations  in  the  early  part  of  the  direct  P  arrivals  at  our 
stations,  which  are  at  upper  mantle  triplication  distances. 
The  P-wave  coda  should  involve  scattered  arrivals,  some 
of  which  originate  in  the  Lg  and  Sn  windows,  so  there 
could  be  some  sensitivity  to  the  path  properties.  Figure 
1 1  shows  the  varying  degree  of  correlation  between  am¬ 
plitudes  of  each  phase  and  the  maximum  water  depth. 
Because  of  the  upper  mantle  distance  ranges  involved, 
no  simple  distance  correction  is  realistic  for  the  P-wave 
windows,  so  in  this  comparison  we  make  no  geometric 
spreading  correction  for  Sn  or  Lg.  Comparison  with  Fig¬ 
ures  9  and  10  indicates  that  the  geometric  spreading  cor¬ 
rections  used  for  the  latter  two  phases  are  actually  not 
very  significant  relative  to  the  trends  found  with  maxi¬ 
mum  water  depth.  Even  without  any  distance  correction, 
Lg  has  the  greatest  slope  and  correlation  coefficient  with 
maximum  depth,  followed  by  Sn ,  and  then  the  0-  to  50- 
sec  P- wave  window,  which  includes  substantial  P  coda. 
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Figure  8.  Empirical  relation  between  the 
rms  amplitude  ratio  of  Sn/Lg  versus  along- 
path  topographic  roughness  (variance)  for 
observations  from  the  Semipalatinsk  source 
region  (right)  and  the  Novaya  Zemlya  source 
region  (left).  The  strong  trend  on  the  right 
was  first  reported  by  Zhang  and  Lay  (1994), 
and  indicates  that  irregularities  in  the  wave 
guide  structure  influence  the  regional  phase 
energy  partitioning.  The  lack  of  a  strong 
trend  on  the  left  indicates  that  additional 
factors,  presumably  associated  with  the  un¬ 
derwater  segments  of  the  paths,  dominate  in 
the  Novaya  Zemlya  data. 
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Figure  9.  Correlation  of  rms  Lg  ampli¬ 
tudes  from  the  Novaya  Zemlya  explosions 
with  four  parameters  that  characterize  the 
underwater  segments  of  each  path,  mean 
depth,  underwater  length,  maximum  depth, 
and  underwater  area.  The  data  are  mean 
values  on  each  path,  after  equalizing  the 
source  strengths  to  a  common  mb  —  5.0,  and 
correcting  for  distance.  The  geometrical 
spreading  factor  for  Lg  is  approximated  by 
the  power  law  r~5/ 6  (Nuttli,  1973). 
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Figure  10.  Similar  to  Figure  9  but  for  Sn 
amplitudes.  The  r"5/6  correction  is  still  used. 
The  results  are  close  to  those  in  Figure  9, 
but  the  correlation  and  trend  with  maximum 
depth  are  a  little  weaker. 
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The  direct  0-  to  5-sec  P- wave  rms  amplitude  has  no  cor¬ 
relation  with  the  water  depth,  as  expected.  Note  that  sta¬ 
tion  TLY,  a  consistent  outlier  for  both  Sn  and  Lg ,  is  also 
low  for  P,  suggesting  a  receiver  site  effect  or  possibly 
an  error  in  the  instrument  corrections  and/or  gains  pro¬ 
vided  in  the  data  exchange. 

Given  the  complex  nature  of  geometric  spreading  and 
path  effects,  we  can  better  explore  energy  partitioning 
by  examing  rms  ratios.  Ratios  are  also  of  importance 
because  ratios  such  as  P/Lg  are  being  explored  as  the 
main  regional  discriminant  between  nuclear  explosions 
and  earthquakes.  We  seek  to  establish  whether  such  ra¬ 
tios  are  affected  by  the  water  segment.  We  expect  that 
ratios  for  regional  phases  may  be  affected  by  average 
path  properties,  as  shown  in  Figure  8,  and  we  can  ex¬ 
plore  whether  applying  corrections  for  the  on-land  por¬ 
tion  of  the  paths  enhances  the  sensitivity  to  the  under¬ 
water  segments.  Figure  12  shows  results  for  the  Sn/Lg 
ratios,  correlated  with  the  four  underwater  parameters. 
The  Sn/Lg  has  a  fairly  strong  relationship  with  maximum 
depth,  which  can  be  attributed  to  the  stronger  trend  in 
Figure  9  versus  Figure  10.  Note  that  the  elimination  of 
any  site  effect  for  TLY  tightens  up  the  trend.  Also  note 
that  the  single  observation  at  FRU  should  be  down¬ 
weighted. 

From  Figure  1  it  is  clear  that  the  on-land  path  seg¬ 
ments  are  much  longer  than  the  underwater  segments. 


We  use  the  correlations  between  Sn/Lg  and  P/Lg  and  on- 
land  topographic  roughness  of  the  path  found  in  our  pre¬ 
vious  study  (Zhang  and  Lay,  1994)  to  make  corrections 
for  on-land  roughness  for  the  Novaya  Zemlya  data.  Since 
the  range  in  path  roughness  is  relatively  small  (see  Fig. 
8),  these  corrections  are  not  very  large,  but  one  should 
observe  an  improvement  in  the  correlations  with  under¬ 
water  properties  if  the  on-land  contributions  are  in  fact 
suppressed.  The  results  for  Sn/Lg  after  correction  for 
roughness  (relative  to  the  path  to  ARU)  are  shown  in 
Figure  13.  The  correlations  of  Sn/Lg  ratios  with  mean 
depth,  maximum  depth,  and  water  area  are  all  actually 
improved,  lending  credibility  to  the  influence  of  the  un¬ 
derwater  segment.  The  strong  effect  associated  with  water 
depth  does  in  fact  appear  to  be  responsible  for  the  poor 
correlation  with  overall  path  roughness  for  these  data 
found  in  Figure  8. 

Since  Lg  shows  a  trend  with  underwater  properties, 
while  the  0-  to  5-sec  P  window  does  not,  the  presence 
of  correlations  in  the  P/Lg  ratio  (Fig.  14)  is  no  surprise. 
Because  the  Semipalatinsk  data  span  comparable  dis¬ 
tance  ranges,  it  is  possible  to  use  the  relationship  found 
by  Zhang  and  Lay  (1994)  for  P/Lg  versus  continental 
path  roughness  to  make  a  correction,  but  this  effectively 
corrects  only  for  the  Lg  wave.  It  is  interesting  to  see  (Fig. 
15)  that  application  of  this  correction  does  in  fact  reduce 
the  scatter  in  the  relationship  with  the  underwater  path 
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Figure  11.  (a)  through  (d)  Comparison 

of  relationships  between  rms  amplitude 
measurements  for  different  phases  in  the 
short-period  signals  and  maximum  water 
depth  on  each  path.  No  distance  corrections 
are  applied  to  any  of  the  phases  in  this  case. 
The  direct  P  window,  measured  for  the  first 
5  sec  of  the  waveform,  is  shown  in  (b),  and 
no  relationship  to  the  water  depth  is  seen. 
A  weak  trend  is  found  for  a  long,  0-  to  50- 
sec  P  window,  which  includes  substantial 
coda,  and  stronger  trends  are  found  for  Sn 
and  Lg. 
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Figure  12.  Correlations  between  Sn/Lg 
ratios  and  the  four  underwater  path  param¬ 
eters. 
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Figure  13.  Similar  to  Figure  12,  but  now 
the  Sn/Lg  ratios  have  been  corrected  for  rel¬ 
ative  path  roughness,  using  the  continental 
path  correlations  found  for  Semipalatinsk 
data  by  Zhang  and  Lay  (1994). 
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Figure  14,  Correlations  between  P/Lg 
and  the  four  underwater  path  parameters. 
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Figure  15.  Similar  to  Figure  14,  but  now 
the  P/Lg  ratios  have  been  corrected  for  rel¬ 
ative  path  roughness,  using  the  continental 
path  correlations  found  for  Semipalatinsk 
data  by  Zhang  and  Lay  (1994). 
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properties,  with  the  best  correlation  again  being  with 
maximum  depth.  Similar  results  are  found  for  the  0-  to 
50-sec  P- wave  window.  While  these  data  are  at  larger 
distances  than  conventionally  used  in  regional  discrim¬ 
ination  efforts  (based  on  Pti/Lg),  these  results  hold  much 
promise  for  comparable  efforts  to  reduce  scatter  in  dis¬ 
criminant  ratios  used  for  nonproliferation  monitoring. 

Discussion  and  Conclusions 

This  study  indicates  that  Lg  energy  is  lost  when  tra¬ 
versing  an  underwater  continental  margin,  possibly  ac¬ 
companied  by  a  corresponding  increase  in  Sn  energy. 
Studies  by  Isacks  and  Stephens  (1975),  Chinn  et  aL 
(1980),  and  Ni  and  Barazangi  (1983)  have  discussed  so- 
called  “early  Lg  phases  in  the  Sn  coda,  which  arrive 
several  seconds  ahead  of  the  3.5  km/sec  group  velocity 
arrival  time.  Baumgardt  (1990)  also  observed  early  Lg 
in  his  study  of  data  at  NORESS  and  Graefenburg  arrays 
for  Novaya  Zemlya  sources.  He  argues  that  the  primary 
mechanism  for  energy  partitioning  is  Lg-to-Sn  and  Sn- to- 
Lg  scattering.  Our  study  emphasizes  the  important  role 
played  by  the  underwater  path,  and  we  agree  that  Lg- to- 
Sn  and  Sn-to-Lg  scattering  probably  plays  some  role  in 
the  propagation.  This  could  explain  the  enhanced  Sn/Lg 
ratios  of  the  Novaya  Zemlya  observations  relative  to 
Semipalatinsk  data.  However,  it  is  important  to  note  that 
the  trends  in  Sn/Lg  and  P/Lg  ratios  versus  bathymetric 
parameters  in  Figures  13  and  15  may  be  primarily  an  Lg 
effect.  Moreover,  it  is  not  possible  for  us  to  confidently 
isolate  the  location  of  any  conversions  that  do  take  place, 
other  than  to  say  that  it  must  be  in  the  general  vicinity 
of  the  source  or  receivers,  given  our  gross  group  velocity 
windowing  procedure. 

The  Barents  Sea  and  Kara  Sea  are  similar  shallow 
basins.  Gramberg  (1988)  identifies  two  defining  char¬ 
acteristics  of  the  central  Barents  basin.  One  is  thinning 
of  the  crust  by  about  10  km  (Fig.  16b).  Another  is  the 
“missing  granitic  layer.”  The  data  that  Baumgardt  (1990) 
analyzes  traverse  the  central  Barents  Sea,  where  signif¬ 
icant  crustal  thinning  is  expected.  In  our  data  set  the  three 
paths  sampling  the  Barents  Sea  cross  the  Pechora  plate, 
in  the  southeast  region  of  the  basin,  where  the  crustal 
thickness  is  about  35  km.  Bogolepov  et  al.  (1990)  stud¬ 
ied  the  structure  of  the  Kara  Sea  basin  with  a  combi¬ 
nation  of  magnetic,  gravitational,  and  seismic  refraction 
data.  They  find  similar  characteristics  in  the  center  of 
south  Kara  syncline  (Fig.  16c)  to  the  central  Barents  Sea. 
The  depth  of  the  Moho  varies  from  32  to  35  km  in  the 
near-shore  region  to  26  to  28  km  in  the  central  Kara 
syncline.  The  depth  of  the  Conrad  discontinuity  varies 
from  30  to  22  to  25  km.  In  the  central  region,  the  crys¬ 
talline  granitic  complex  is  replaced  by  a  gabbroic  com¬ 
plex.  Another  characteristic  is  the  complex  geometry  of 
the  granitic  basement,  with  extensive  folding,  meta¬ 
morphism  and  imbricate  structure.  Four  of  our  paths,  to 


BOD,  NRI,  TLY,  and  NVS  cross  the  central  Kara  Sea 
basin,  while  the  last  two,  to  FRU  and  ARU,  traverse  the 
Ural-Novaya  Zemlya  fold  belt.  We  observe  significant 
reductions  in  Lg  amplitudes  for  the  paths  across  the  Kara 
Sea. 

The  Lg  phase  can  be  viewed  as  a  guided  wave  in  the 
crust,  very  dependent  upon  the  lateral  continuity  of  the 
wave  guide.  Given  that  disruptions  of  the  wave  guide 
are  expected  near  continental  margins,  and  even  within 
tectonically  active  continental  regions,  there  are  many 
concerns  about  the  reliability  of  Lg  in  discrimination 
measurements  (e.g.,  Kennett,  1989),  independent  of  its 
remarkable  relative  stability  on  a  given  path  for  relative 
yield  estimation  (e.g.,  Ringdal  et  al.y  1992;  Hansen  et 
al .,  1990).  Crustal  thinning  apparently  can  reduce  Lg 
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Figure  16.  (a)  Map  of  the  adjacent  region  of 

the  islands  of  Novaya  Zemlya,  showing  two  lines 
of  cross  section  and  eight  paths,  (b)  NW-SE  cross 
section,  labeled  A-B  in  (a),  across  the  southeast¬ 
ern  part  of  the  Barents  Sea  basin  (after  Gramberg, 
1988).  (c)  NE-SW  cross  section  of  the  Kara  Sea 
basin  (after  Bogolepov  et  al .,  1990),  labeled  N-S 
in  (a),  across  the  North  Kara  plate  and  south  Kara 
syncline. 
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amplitude  (Cao  and  Muirhead,  1993).  Thus,  we  may  ex¬ 
plain  the  systematic  trends  between  the  Barents  group 
and  Kara  group  as  a  result  of  variable  degree  of  blockage 
associated  with  thinning  of  the  crustal  wave  guide.  How¬ 
ever,  the  precise  mechanisms  involved  are  subject  to  de¬ 
bate. 

It  is  quite  possible  that  sediment  thickness  variations 
in  the  basins  influence  Lg  blockage,  but  we  do  not  have 
a  robust  procedure  for  separating  the  effects  of  crustal 
thinning  and  sediment  thickness,  which  are  strongly  cor¬ 
related.  Maximum  water  depth,  the  parameter  yielding 
the  strongest  relationships  with  Lg  amplitudes,  is  strongly 
affected  by  near-surface  geology,  such  as  the  sediment 
thickness  in  the  basins.  Some  information  is  available 
about  the  thickness  of  sedimentary  basins,  along  with 
Moho  depth  along  the  paths  sampled  by  our  data,  mainly 
from  Deep  Seismic  Sounding  and  seismic  reflection  pro¬ 
filing  as  well  as  other  geological  and  geophysical  data 
(Kunin  and  Sheykh-Zade,  1983;  Kunin,  1986).  This  in¬ 
formation  has  been  digitized  and  made  available  by  the 
Cornell  Institute  for  Studies  of  the  Continents,  with  0.25° 
resolution.  Figure  17  shows  the  relationships  between  Lg 
amplitudes,  corrected  for  geometric  spreading,  and  max¬ 
imum  water  depth  (Fig.  17a),  minimum  crustal  thickness 
(Fig.  17b),  and  maximum  sedimentary  layer  thickness 
(Fig.  17c)  along  each  path.  While  the  data  set  is  limited, 
the  strongest  variation  is  found  with  water  depth,  and  no 


trend  is  apparent  with  minimum  crustal  thickness,  which 
basically  measures  the  necking  of  the  crust  along  the 
continental  margin.  However,  the  latter  parameter  is  not 
well  sampled  by  our  particular  data  set.  For  other  data 
sets  (Zhang,  Schwartz,  and  Lay,  1994),  crustal  necking 
does  appear  to  strongly  affect  Lg  amplitudes.  Sediment 
layer  thickness  is  associated  with  a  weak  trend,  but  the 
statistical  significance  is  low.  The  precise  mechanism  by 
which  the  Lg  energy  is  lost  requires  further  study.  The 
complex  structure  of  the  crystalline  basement  may  ac¬ 
count  for  some  of  the  observed  behavior,  such  as  the  low 
amplitudes  of  both  Lg  and  Sn  along  the  path  to  TLY  (Fig. 
10- 

The  foregoing  discussion  has  emphasized  the  pos¬ 
sible  elastic  scattering  effects  of  the  disrupted  wave  guide. 
There  is  no  question  that  attenuation  is  very  important 
for  Lg  and  Sn  phases  as  well,  and  it  would  be  very  at¬ 
tractive  to  correct  our  data  for  attenuation  variations  on 
the  different  paths.  Many  efforts  have  evaluated  fre¬ 
quency-dependent  attenuation  of  Lg  in  western  Europe 
(e.g.,  Nicolas  et  al. ,  1982;  Campillo  et  al.,  1985;  Cam- 
pillo,  1987)  and  Eurasia  (Nuttli,  1980,  1981;  Chun  et 
al .,  1992).  Analysis  of  Lg  coda  also  provides  constraints 
on  the  Lg  quality  factor  (e.g.,  Xie  and  Mitchell,  1990, 
1991;  Pan  et  al. ,  1992;  Mitchell  et  al. ,  1993).  The  Lg 
coda  Q  values  agree  well  with  Lg  Q  values  in  some  cases 
(Mitchell  et  al. ,  1993).  We  use  the  Eurasian  Lg  coda  Q 


Figure  17.  The  relationships  between  Lg 
amplitudes,  corrected  for  geometric  spread¬ 
ing,  and  (a)  maximum  water  depth,  (b) 
minimum  crustal  thickness,  (c)  maximum 
sedimentary  layer  thickness,  and  (d)  the  at¬ 
tenuation  term  y\  along  each  path. 
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map  of  Pan  et  al.  (1992)  to  determine  the  exponential 
attenuation  term  yA  (Nuttli,  1986)  for  a  frequency  of  1 
Hz  for  each  of  our  paths  (Fig.  17d).  There  is  a  clear 
correlation  with  log  Lgi  as  expected,  but  the  trend  is  quite 
a  bit  weaker  than  expected  (the  expected  slope  is  -log 
(e)  =  -0.43).  If  we  apply  the  attenuation  corrections, 
using  the  simplifying  approximation  that  our  rms  values 
are  representative  of  the  1-Hz  amplitudes,  the  Lg  ampli¬ 
tudes  no  longer  show  any  trend  with  bathymetric  prop¬ 
erties.  However,  our  feeling  is  that  the  Q  model  has  the 
right  geographic  pattern,  but  it  overpredicts  attenuation 
variations  on  the  different  paths,  when  applied  to  Lg  di¬ 
rectly,  as  suggested  by  the  low  slope  seen  in  Figure  17d. 
Thus,  we  have  not  used  corrected  Lg  amplitudes  in  our 
attempt  to  detect  underwater  path  effects.  We  found  that 
Sn  amplitudes  show  a  strong  trend  with  the  Lg  coda  Q 
variations  for  Novaya  Zemlya,  but  the  Sg/Lg  ratios  do 
not  vary  with  yA  for  either  test  site,  so  it  is  not  likely 
that  attenuation  plays  a  key  role  in  the  strong  patterns 
seen  in  Figure  8.  It  is  not  clear  how  to  correct  the  Sn 
data  for  attenuation,  but  there  is  a  clear  suggestion  that 
Sn  attenuation  is  linked  to  the  Lg  coda  Q  variations.  This 
suggests  coupling  of  upper  mantle  and  crustal  attenua¬ 
tion,  and  raises  the  possibility  that  the  Lg  coda  Q  model 
may  include  a  mantle  component  that  overcorrects  the 
Lg  phases.  Until  we  have  independent  models  for  Sn  at¬ 
tenuation  we  cannot  proceed  further,  but  it  is  clear  that 
attenuation  is  an  important  factor  influencing  the  re¬ 
gional  phase-amplitude  ratios,  and  significant  progress 
is  being  made  in  mapping  out  the  attenuation  variations 
(Pan  et  al .,  1992). 

Since  we  generally  lack  detailed  knowledge  of  crus¬ 
tal  structure  for  paths  traversed  by  regional  phases  used 
in  nonproliferation  monitoring,  we  have  explored  the 
possibility  of  developing  empirical  path  calibration  pro¬ 
cedures.  Water  depth  proved  to  be  the  most  effective 
elastic  path  characteristic  for  the  Novaya  Zemlya  data 
set,  presumably  since  it  serves  as  a  surrogate  for  the  deep 
structure  of  the  shallow  sea  basin.  Relative  to  the  ref¬ 
erence  data  set  provided  by  Semipalatinsk  recordings, 
almost  all  Novaya  Zemlya  data  exhibit  a  shift  of  the  peak 
energy  toward  the  Sn  window,  whether  the  maximum 
water  depth  is  100  or  300  m.  This  indicates  that  an  am¬ 
plitude  correction  is  generally  needed  for  wave  paths  that 
cross  marginal  seas.  In  addition,  200-m  water-depth  dif¬ 
ferences  can  reflect  variations  in  crustal  thickness  and 
sediment  thickness,  causing  Lg  amplitudes  to  vary.  Us¬ 
ing  empirical  corrections  should  enable  significant  re¬ 
duction  in  the  scatter  of  both  Sn/Lg  and  P/Lg  amplitude 
ratios.  Combined  with  empirical  corrections  for  topog¬ 
raphy  along  the  on-land  paths  and  corrections  for  lat¬ 
erally  varying  attenuation  models  in  the  crust  and  upper 
mantle,  reduced  scatter  in  such  measurements  may  im¬ 
prove  nonproliferation  monitoring  efforts  in  diverse  crustal 
environments. 
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Abstract  Four  data  sets  characterizing  gross  crustal  waveguide  configuration  and  attenuation 
properties  in  Eurasia  (surface  topography,  Moho  depth,  sediment  thickness,  and  L  coda  Q 
value)  are  used  to  examine  path  influences  on  short-period  regional  P/Lg  ratios.  Linear 
regressions  show  considerable  correlations  between  log  P  lLg  and  waveguide  properties  for 
both  earthquake  and  explosion  data.  This  is  of  interest  because  P/Lg  ratios  are  considered  to 
be  promising  regional  discriminants  for  which  it  is  desirable  to  reduce  scatter  caused  by  propa¬ 
gation  effects.  To  develop  a  comprehensive  and  stable  model  describing  the  path  effects,  mul¬ 
tivariate  regression  is  applied  to  the  data  set.  The  waveguide  properties  considered  are  highly 
correlated  with  each  other,  causing  collinearity  in  the  regression  process.  After  backward 
elimination,  we  obtain  a  final  model  with  statistically  significant  regression  coefficients,  which 
involves  two  independent  variables:  the  path  attenuation  term  yA  and  maximum  sediment 
thickness  on  each  path.  These  two  factors  are  associated  with  (1)  the  overall  intrinsic  and 
scattering  attenuation  in  the  waveguide  and  (2)  the  localized  blockage  effects  caused  by 
waveguide  geometry,  respectively.  Using  the  final  model  to  correct  the  observed  P II  meas¬ 
urements  reduces  the  variances  for  separate  earthquake  data  and  explosion  data  by  40%  and 
27%,  respectively.  This  correction  slightly  enhances  the  performance  of  the  discriminant  for 
periods  near  Is. 


Introduction 

The  manner  in  which  regional  Pn ,  Pgi  Sny  and  Lg 
phases  are  affected  by  passing  through  different  geologi¬ 
cal  structures  is  not  clearly  established  in  terms  of  the 
geometrical  characteristics  and  internal  properties  of  the 
crust.  Since  Lg  is  a  guided  wave  in  the  crustal 
waveguide,  lateral  variations  of  the  waveguide  structure 
should  affect  its  propagation.  Either  thinning  or  thicken¬ 
ing  of  the  waveguide  should  result  in  a  portion  of  the  Lg 
mode  energy  leaking  out,  as  illustrated  using  ray 
diagrams  by  Kennett  [1986].  However,  it  has  been  hard 
to  verify  this  qualitative  inference  with  numerical  simu¬ 
lations.  Campiilo  [1987]  compared  synthetic  seismo¬ 
grams  computed  in  crustal  models  with  variations  of 
Moho  depth  and  sediment  thickness  with  those  computed 
in  a  flat-layered  model.  He  found  that  the  wave  shapes 
are  sensitive  to  the  presence  of  the  irregularities  but  the 
overall  amplitude  decay  was  not.  Maupin  [1989]  per¬ 
formed  numerical  modeling  of  Lg  propagation  in  a 
simplified  model  of  the  North  Sea  graben.  This  model 
did  not  predict  the  severe  Lg  attenuation  observed  in  this 
area.  On  the  contrary,  the  synthesized  Lg  wave  train 
appears  surprisingly  robust  when  crossing  a  zone  where 
its  waveguide  is  strongly  deformed.  Regan  and  Har- 
krider  [1989]  simulated  the  extreme  cases  of  transitions 
from  continent  to  ocean  crust  and  vice  versa,  for 
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transverse-component  Lg  propagation.  In  this  case,  the 
Lg  amplitude  did  decrease,  but  the  attenuation  is  not 
sufficient  to  explain  observed  values.  The  strongest 
attenuation  predicted  involves  a  reduction  of  a  factor  of 
2-3  in  amplitude.  However,  it  is  well  known  that  Lg  is 
totally  eliminated  when  traversing  oceanic  crust  [Press 
and  Ewing ,  1952].  Thus  other  factors  and/or  more  com¬ 
plex  structures  must  be  considered  to  explain  the 
observed  attenuation  of  Lg.  Recently,  Cao  and  Muir- 
head  [1993]  successfully  modeled  Lg  blockage  with  a 
P—SV  finite  difference  method  after  introducing  a  water 
layer  above  the  free  surface.  The  water  layer  strongly 
attenuates  Rg1  the  short-period  fundamental  mode  Ray¬ 
leigh  wave,  but  has  less  effect  on  Lg.  However, 
interpretation  of  the  influence  of  variable  waveguide 
structure  without  water-covered  segments  is  still  open. 

The  most  important  factor  with  predictable  effect  on 
Lg  other  than  the  waveguide  structure  is  the  quality  fac¬ 
tor  Q,  accounting  for  both  intrinsic  attenuation  and 
scattering  losses.  Many  investigations  of  Q  using  the 
decay  of  Lg  amplitude  with  epicentral  distance  have 
been  conducted  in  western  Europe  [e.g.,  Nicolas  et  al. , 
1982;  Campiilo  et  al. ,  1985]  and  Eurasia  [e.g.,  Nuttli , 
1980,  1981;  Chun  et  al. ,  1992].  Another  way  to  evalu¬ 
ate  waveguide  Q  is  from  frequency  dependent  decay  of 
Lg  coda  [e.g.,  Xie  and  Mitchell ,  1990;  Pan  et  al. ,  1992; 
Mitchell  et  al ,  1993].  The  Q  obtained  in  this  way  is 
called  Lg  coda  Q .  Lg  Q  values  sometimes  agree  well 
with  Lg  coda  Q  values  in  the  same  region  [Mitchell  and 
Hwang ,  1987;  Mitchell  et  al ,  1993].  However,  they 
may  have  significantly  different  values  in  some  regions 
[Campiilo,  1990].  Both  geometrical  blockage  and 
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waveguide  attenuation  affect  Lg,  but  it  is  difficult  to 
clearly  separate  these  effects  observationally.  Indeed,  dis¬ 
tortions  of  the  waveguide  are  likely  to  be  associated  with 
tectonic  factors  influencing  intrinsic  attenuation,  so  the 
effects  may  often  be  coupled. 

It  is' desirable  to  quantitatively  identify  the  respective 
roles  played  by  large-scale  structural  heterogeneities  and 
by  the  attenuation  properties  of  the  crust,  but  at  present, 
doing  this  requires  empirical  approaches.  Because  of  the 
lack  of  information  about  the  lower  crust  and  upper 
mantle,  Zhang  and  Lay  [1994a]  first  used  surface  topog¬ 
raphy  as  a  manifestation  of  the  varying  crustal  structure. 
They  found  a  surprisingly  strong  correlation  between 
SnILg  ratios  for  Eurasian  explosions  and  roughness  or 
mean  altitude  of  the  topography  along  the  path  to  each 
station.  Zhang  and  Lay  [1994b]  found  that  the  maximum 
water  depth  along  each  path  correlates  with  the  Lg  and 
Sn  amplitudes  for  signals  traversing  the  Kara  and 
Barents  Seas.  Jih  [1993]  has  demonstrated  strong 
scattering  of  Rg  phases  by  surface  topography  which 
explains  why  Rg  seldom  propagates  to  large  distances. 
Rg  to  Lg  scattering  may  produce  some  correlation 
between  Lg  behavior  and  surface  topography  characteris¬ 
tics. 

Topography  only  indirectly  reflects  the  properties  of 
the  crust;  to  further  explore  the  nature  of  regional 
seismic  wave  propagation  in  crustal  waveguides  in  this 
paper,  we  use  four  data  sets  that  grossly  characterize  the 
crustal  structure  and  attenuation  properties:  surface 
topography,  Moho  depth,  sediment  thickness,  and  Lg 
coda  Q  values.  Although  each  characterization  of  the 
crust  has  limited  resolution,  complete  representations  are 
available  for  Eurasia.  We  seek  a  model  to  define  the 
relative  importance  of  various  waveguide  characteristics 
for  regional  phases  using  P!Lg  amplitude  ratios  for 
earthquake  and  explosion  recordings  in  Eurasia.  Only  a 
few  such  data  have  overlapping  paths  in  Eurasia,  so  we 
treat  large  populations  of  paths  rather  than  focusing  on  a 
single  region.  Not  surprisingly,  many  factors  appear  to 
influence  the  data,  so  multivariate  analysis  is  used  to 
assess  their  importance  and  correlation.  For  example,  in 
a  study  of  regional  wave  propagation  in  northern 
Eurasia,  weak  Lg  amplitude  is  related  to  thinner  crust 
[Zhang  and  Lay ,  1994b].  However,  in  the  Tibetan  Pla¬ 
teau,  where  the  Lg  amplitudes  are  also  weak,  the  crust  is 
the  thickest  in  the  world.  If  the  low  crustal  Q  in  the 
Tibetan  Plateau  is  accounted  for,  the  low  Lg  amplitudes 
can  be  explained.  A  multivariate  description  is  thus 
needed  to  predict  regional  phase  effects.  A  first  genera¬ 
tion  model  is  obtained  here. 

P!Lg  amplitude  or  spectral  ratios  are  the  most  promis¬ 
ing  regional  distance  seismic  discriminants.  However, 
strong  path  effects  on  the  ratios  have  so  far  restricted 
their  use  to  nearly  common  propagation  paths 
[Baumgardt  and  Young ,  1990;  Bennett  et  al ,  1992]. 
Reduction  in  P/Lg  variance  by  correcting  this  ratio  for 
empirically  derived  propagation  effects  may  remove  the 
restriction  of  common  paths  for  explosion  and  earth¬ 
quake  data.  This  would  greatly  enhance  the  possibility  of 
discrimination  in  the  context  of  nonproliferation  monitor¬ 
ing,  where  no  prior  explosion  population  may  be  avail¬ 
able.  This  provides  further  motivation  for  the  empirical 
approach  of  this  paper. 


Data 

Our  explosion  data  are  recordings  from  Soviet-run  sta¬ 
tions  that  were  collected  and  digitized  as  part  of  a  data 
exchange  associated  with  the  Joint  Verification  Experi¬ 
ment.  We  use  211  recordings  of  77  underground  nuclear 
tests  located  in  the  Semipalatinsk  test  site,  a  subset  of 
the  data  used  by  Zhang  and  Lay  [1994a]  after  exclusion 
of  the  data  lacking  P  or  Lg.  Seventy-two  recordings  of 
20  underground  nuclear  tests  located  in  the  Novaya 
Zemlya  test  site  near  Matochkin  Shar  are  extracted  from 
the  data  used  by  Zhang  and  Lay  [1994b].  The  earthquake 
data  set  consists  of  34  digital  waveforms  for  17  events 
between  1988  and  1992  recorded  at  IRIS/IDA  (Incor¬ 
porated  Research  Institutions  for  Seismology/Interna¬ 
tional  Deployment  of  Accelerometers),  CDSN  (Chinese 
Digital  Seismograph  Networks)  and  ASRO  (Abbreviated 
Seismic  Research  Observatory)  stations  in  Eurasia. 
Almost  all  of  the  explosion  data  are  in  the  distance  range 
14°  to  36°,  while  the  earthquake  data  range  from  4°  to 
35°.  The  explosion  data  are  short-period,  while  the  earth¬ 
quake  data  are  mostly  broadband.  For  all  data,  the  instru¬ 
ment  responses  are  equalized  to  the  1988  response  of  the 
CKM-3  seismometer  at  station  OBN,  and  waveforms  are 
bandpass  filtered  from  0.6  to  3.0  Hz  (following  the  pro¬ 
cedure  of  Israelsson  [1992]).  RMS  P  wave  amplitudes 
are  measured  for  a  0-50  s  time  window  (following  Ben - 
nett  et  al.  [1992]),  and  RMS  Lg  amplitudes  are  calcu¬ 
lated  for  a  3.7-3. 1  km/s  group  velocity  window  [Israels¬ 
son,  1992],  A  noise  correction  is  applied  to  each  RMS 
measurement  [Zhang  and  Lay ,  1994a].  We  average  the 
log  P/Lg  ratio  values  of  the  explosion  data  along  each 
path,  yielding  16  path-averaged  explosion  log  P!Lg 
values. 

The  surface  topography  data  are  from  the  5  arc  min 
interval  ETOP05  data  set  compiled  by  the  National  Geo¬ 
physical  Data  Center,  Boulder,  Colorado.  Figure  1  shows 
our  explosion  data  path  coverage  superimposed  on  the 
topography  of  northern  Eurasia.  The  earthquake  data 
path  coverage  is  shown  in  Figure  2  superimposed  on  the 
crustal  thickness  information.  Unfortunately,  the  explo¬ 
sion  and  earthquake  data  do  not  sample  the  same  region. 
The  Moho  depth  (Figure  2)  and  sediment  thickness  (Fig¬ 
ure  3)  information  were  compiled  by  Kunin  and  Sheykh - 
Zade  [1983],  and  Kunin  [1987]  from  a  large  number  of 
deep  seismic  sounding  (DSS)  profiles  and  other  data.  The 
compiled  crustal  parameters  were  digitized  and  gridded 
at  Cornell  University  using  a  1/4  degree  interval  [Field¬ 
ing  et  al. ,  1992].  The  accuracy  of  these  maps  is  difficult 
to  assess,  but  at  least  they  provide  a  guide  as  to  gross 
crustal  structure. 

Lg  coda  Q  values  across  Eurasia  (Figure  4)  were  com¬ 
puted  by  tomographic  inversion  of  explosion  and  earth¬ 
quake  data  by  Pan  et  al.  [1992].  We  use  this  model  as 
a  guide  to  gross  attenuative  properties  of  the  crust, 
recognizing  that  some  effects  of  waveguide  irregularity 
may  be  folded  into  the  Q  values. 

For  each  path  from  source  to  receiver,  great  circle 
variations  of  each  of  the  four  crustal  properties  are 
analyzed.  Examples  of  variations  in  these  profiles  along 
three  different  paths  from  an  earthquake  (September  17, 
1989)  in  the  middle  of  the  Caspian  Sea  to  stations  OBN, 
ARU,  and  ANTO  are  shown  in  Figure  5.  Topography, 
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Figure  1.  Map  showing  the  locations  of  the  former  Soviet  Union  test  sites  at  Novaya  Zemlya  and 
Semipalatinsk,  marked  with  triangles,  and  the  seismological  stations,  marked  with  circles  and  indi¬ 
cated  by  their  codes.  The  lines  between  the  triangles  and  the  circles  correspond  to  great  circle  paths. 
Topography  is  indicated  with  contour  intervals  every  1000  m  from  -2000  to  3000  m  above  sea  level. 
In  order  to  outline  the  Ural  Mountains  and  other  low  plateaus,  500-m  contours  are  also  shown.  Map 
software  is  from  Wessel  and  Smith  [1991]. 


Moho  depth,  and  sediment  thickness  are  drawn  together 
on  profiles  with  a  10-km  lateral  sampling  interval  for 
each  path  (Figure  5a).  The  Lg  coda  Q  profiles  are  shown 
in  Figure  5b  with  a  100-km  lateral  sampling  interval. 
Since  deterministic  modeling  of  regional  phases  along 
these  complex  waveguides  is  not  currently  viable,  we 
parameterize  the  path  properties  and  attempt  to  relate 
summary  parameters  to  the  observed  wave  field  charac¬ 
teristics,  as  manifested  in  the  logarithmic  ratios  of  RMS 
P/RMS  Lg . 

Linear  Regression  of  Waveguide  Parameters 

We  first  explore  the  individual  correlations  between 
waveguide  parameters  and  observed  PILg  amplitude 
ratios  for  our  explosion  and  earthquake  data  sets.  The 
results  for  earthquake  data  are  summarized  in  Figures  6 
and  7.  Those  for  Novaya  Zemlya  explosion  data  are  in 
Figures  8  and  9.  Those  for  Semipalatinsk  explosion  data 
are  in  Figures  10  and  1 1 . 

Topography 

Surface  topography  variations  are  small  relative  to 
crustal  thickness  but  to  some  degree  do  reflect  lluttua- 


tions  in  crustal  thickness  through  isostasy.  Topography 
is  also  affected  by  lateral  variations  in  thermal  structure 
of  the  crust  and  upper  mantle.  These  considerations  sug¬ 
gest  that  topography  may  influence  regional  phases  both 
directly  as  a  source  of  scattering  and  indirectly  as  a  man¬ 
ifestation  of  waveguide  structure  and  thermal  variations. 
Surface  topography  is  also  widely  available  with  high 
precision,  so  it  is  important  to  establish  its  influence,  if 
any,  on  regional  phases. 

We  examine  the  correlation  of  log  P/Lg  ratio  along 
each  source-receiver  path  with  corresponding  values  of 
both  minimum  and  mean  altitude.  The  results  for  earth¬ 
quakes  (Figures  6a  and  6b),  and  for  explosions  at 
Novaya  Zemlya  (Figures  8a  and  8b).  and  Semipalatinsk 
(Figures  10a  and  10b)  all  indicate  similar  trends  with 
decreasing  log  P/Lg  ratios  as  minimum  or  mean  altitude 
increases.  All  correlations  are  statistically  significant 
with  the  exception  of  mean  altitude  for  the  Novaya  Zem¬ 
lya  data,  which  has  the  smallest  range  in  mean  altitude 
among  the  three  data  sets.  Because  of  the  smaller  range 
in  topography  sampled  by  both  explosion  data  sets  com¬ 
pared  with  the  earthquakes,  the  slope  of  the  best  fit  line 
is  most  reliably  determined  by  the  earthquakes. 

The  results  in  Figure  6  are  the  first  to  indicate  statisti¬ 
cal  topographic  influence  on  regional  phases  from  earth- 
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Figure  2.  Map  showing  the  locations  of  17  earthquakes  in  Central  Asia,  the  Caucasus,  and  the  Mid¬ 
dle  East,  marked  with  triangles,  and  the  seismological  stations  used,  marked  with  circles.  Contours 
indicate  variations  in  Moho  depth,  digitized  at  Cornell,  from  a  Soviet  map  compiled  by  N.  Kunin 
[Kunin  and  Sheykh-Zade ,  1983;  Fielding  et  ai  ,  1992].  Moho  topography  is  indicated  with  10-km 
contour  intervals. 


quakes.  This  pattern  is  consistent  with  the  previous 
results  found  by  Zhang  and  Lay  [1994b]  using  explosion 
data  from  Novaya  Zemlya,  although  they  measured  the 
RMS  P  wave  amplitude  in  a  0-5  s  window  rather  than 
our  0-50  s  window  and  they  correlated  this  with  max¬ 
imum  water  depth  rather  than  mean  and  minimum  alti¬ 
tude.  The  observed  trends  are  most  reasonably  attributed 
to  Lg  propagation,  since  much  of  the  P  wave  energy  at 
these  distances  dives  into  the  upper  mantle.  Zhang  and 
Lay  [1994a]  considered  topographic  influences  on 
regional  phases  in  the  Semipalatinsk  data,  finding  results 
like  those  in  Figure  10b,  although  we  include  an  addi¬ 
tional  station,  ARU,  here  which  does  reduce  the  correla¬ 
tion.  The  earlier  study  limited  the  distance  range  used 
whereas  here  we  include  all  observations.  The  Zhang 
and  Lay  [1994a]  study  reveals  strong  trends  for  SnIL s 
ratios  versus  mean  altitude,  as  well  as  for  measures  of 
average  roughness  of  the  path.  Zhang  and  Lay  [1994b] 
indicated  that  the  variability  associated  with  the 
differences  in  underwater  segments  appears  to 
overwhelm  any  effects  of  surface  roughness  for  Novaya 
Zemlya  data.  In  this  study  we  find  that  the  correlations 
of  log  P/Lg  with  the  roughness  of  surface  topography 
and  variance  of  crustal  thickness  are  also  poor  for  the 
earthquake  data  (not  shown).  This  is  somewhat  surprising 
as  one  might  expect  blockage  effects  to  be  controlled  by 
variability  of  the  crustal  thickness.  The  reason  will  be 


explored  in  future  studies,  along  with  alternate  parame- 
terizations  of  crustal  characteristics. 

Crustal  Thickness 

Crustal  thickness  is  an  important  property  of  the  cru¬ 
stal  waveguide  expected  to  profoundly  affect  propagation 
of  regional  waves.  Thinning  of  the  crust  can  strongly 
reduce  Lg  amplitude  [Cao  and  Muirhead ,  1993;  Zhang 
and  Lay ,  1994b].  We  correlate  both  minimum  and  mean 
crustal  thickness  with  log  P  !L^  ratios  to  explore  whether 
variations  in  L g  amplitude  are  due  to  average  crustal 
structure  or  maximum  thinning  at  restricted  locations. 
The  crust  thins  dramatically  in  several  regions  of  our 
study  area.  In  particular,  the  paths  from  Novaya  Zemlya 
to  Eurasian  stations  traverse  the  thin  crust  of  the  Barents 
and  Kara  Seas  (Figure  1).  The  thinnest  part  of  the  crust 
below  the  Barents  Sea  occurs  near  50°E,  70°N,  where 
the  thickness  is  less  than  30  km.  To  the  south,  the  great 
circle  paths  connecting  earthquakes  to  stations  cross  four 
thin  crustal  regions,  the  southern,  central,  and  northern 
Caspian  Seas  and  the  Aral  Sea  (Figure  2). 

Correlations  of  log  P !LK  with  crustal  thickness  for 
earthquake  data  (Figures  6c  and  6d),  and  for  explosions 
from  Novaya  Zemlya  (Figures  8c  and  8d),  and  Semipala¬ 
tinsk  (Figures  l()c  and  lOd)  are  generally  stronger  than 
those  with  topography,  suggesting  that  direct  measures  of 
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Figure  4.  Map  showing  the  Lg  coda  Q  values  of  Eurasia  at  /  =  1  Hz.  These  tomographic  results 
are  from  Pan  et  al  [1992]. 
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Figure  5.  Examples  of  great  circle  characteristics  of  the 
path  variables  to  three  stations  for  event  September  17, 
1989,  which  occurred  in  the  middle  of  the  Caspian  Sea, 
(a)  The  three  curves  in  each  panel  indicate  altitude,  sedi¬ 
ment  depth,  and  Moho  depth  from  top  to  bottom  at  a 
10-km  horizontal  interval,  (b)  The  curves  indicate  Lg 
coda  Q  values  at  a  100-km  horizontal  interval. 


crustal  thickness  better  reflect  the  nature  of  the 
waveguide  than  surface  topography.  The  shortcoming  of 
crustal  thickness  variables,  however,  is  that  their  accu¬ 
racy  and  reliability  are  not  as  good  as  for  topography, 
especially  in  remote  areas  or  beneath  the  seas. 

Sediment  Thickness 

Sediment  thickness  variations  are  believed  to  have 
strong  blockage  effects  on  Lg  and  perhaps  Sn  propaga¬ 
tion  [Mitchell  and  Hwang,  1987;  Baumgardt ,  1990]. 
Our  study  area  samples  a  few  thick  sedimentary  basins 
(Figure  3),  such  as  in  the  South  Caspian  depression  (sed¬ 
iment  thickness  of  25-30  km)  and  South  Barents  depres¬ 


sion  (20-23  km).  Such  large  sedimentary  sections  have 
not  been  identified  in  other  continents  or  beneath  the 
oceans  [ Kunin ,  1987].  On  the  other  hand,  much  of  the 
remaining  area  in  Eurasia  has  little  to  no  sediment  cover. 
Thus  we  might  expect  large  variations  in  path  effect.  As 
for  crustal  thickness,  the  accuracy  and  reliability  of  sedi¬ 
ment  thickness  data  are  not  as  good  as  for  surface  topog¬ 
raphy. 

We  find  relatively  strong  correlations  between 
log  PILg  ratios  and  along  path  sediment  thickness  for 
earthquake  (Figures  7a  and  7b)  and  explosion  data  (Fig¬ 
ures  9a,  9b,  I  la,  and  lib).  We  evaluated  both  max¬ 
imum  and  mean  sediment  thickness  and  both  show  com¬ 
parable  correlations.  The  slopes  are  consistent  with 
reduced  Lg  amplitudes  for  thicker  sediments  along  the 
path. 


Attenuation 


The  Lg  coda  Q  distributions  for  the  northern  and 
southern  regions  of  the  study  area  are  quite  different. 
High  <2  values  are  found  in  the  north  and  much  lower 
values  in  the  south  (see  Figure  4).  The  lowest  value  is 
found  south  of  the  Caspian  Sea.  The  regional  differences 
make  the  attenuative  properties  for  our  three  data  sets 
very  different.  Without  correction  for  variations  in 
attenuation,  it  is  unlikely  that  discrimination  can  be 
achieved  in  a  region  with  a  wide  range  of  crustal  Q . 

Figures  7c,  7d,  9c,  9d,  11c,  and  lid  show  correlations 
of  log  P ILg  ratio  with  minimum  Q  value  and  the  path- 
integral  yA  value.  The  definition  of  y  follows  Nuttli 
[1986]: 


7 (/) 


UQ(f)' 


where  U  is  the  group  velocity,  taken  as  3.5  km/s,  /  is 
the  wave  frequency,  taken  as  1  Hz,  and 
Q(f)  =  Qo/11-  We  use  the  y  value  at  /  =  1  Hz;  thus 
the  Lg  coda  Q  at  I -Hz  frequency  given  by  Pan  et  al 
[1992]  is  used  directly.  A  is  epicentral  distance  in 
kilometers. 

The  yA  value  has  a  strong  correlation  with  log  P t Lg 
ratio  for  earthquake  data,  but  the  correlations  are  not  as 
significant  for  explosion  data.  In  part  this  may  again  be 
due  to  the  reduced  range  of  yA  sampled  by  the  explosion 
paths.  Semipalatinsk  data  cluster  near  yA  =  4,  which 
may  preclude  finding  any  correlation.  The  minimum  Q 
variable  (Figures  7c,  9c,  and  11c)  also  has  moderate 
correlations  with  log  PILg  ratios.  Localized  regions  of 
strong  attenuation  may  effectively  reduce  Lg  amplitudes, 
but  the  earthquake  data  suggest  that  overall  path  attenua¬ 
tion  has  stronger  correlation  with  log  P/Lg  ratios. 

These  linear  regressions  indicate  that  all  four  crustal 
characterizations  are  at  least  weakly  correlated  with  the 
log  PILg  ratio.  We  would  like  to  develop  path  correc¬ 
tions  to  reduce  the  variance  in  the  measured  ratios,  as 
this  may  be  key  to  effective  discrimination.  However, 
we  need  to  consider  which  parameters  to  use  in  develop¬ 
ing  any  path  corrections.  The  four  crustal  factors  con¬ 
sidered  have  rather  strong  correlations  amongst  them¬ 
selves.  It  is  possible  that  a  limited  number  of  factors 
play  the  main  role  in  the  path  effect,  with  others  corre¬ 
lating  with  log  P/Lg  mathematically,  but  not  physically. 
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Figure  6.  The  log  P!Lg  ratio  variations  with  (a)  minimum  and  (b)  mean  altitude,  and  (c)  minimum 
and  (d)  mean  crustal  thickness  for  earthquake  data.  CC  stands  for  correlation  coefficient,  SIG  is 
standard  deviation  of  linear  correlation,  SLO  is  the  slope  of  the  linear  regression.  "P0-50"  refers  to 
RMS  amplitude  in  50-s  window  beginning  with  P  onset. 


In  order  to  evaluate  each  factor,  to  obtain  a  reasonable 
and  stable  model  for  the  overall  path  effect,  we  use  mul¬ 
tivariate  regression. 


Multivariate  Analysis 

We  have  shown  that  topography,  crustal  thickness, 
sediment  thickness,  and  yA  computed  from  Lg  coda  Q 
and  distance  have  correlations  with  log  P/Lg  ratios.  A 
multivariate  model  is  now  developed  for  the  relationship 
between  the  log  P/Lg  ratios  and  these  factors.  For  topog¬ 
raphy  and  crustal  thickness,  we  include  both  minimum 
and  mean  values  for  the  purpose  of  identifying  whether 
local  properties  or  mean  properties  play  a  more  impor¬ 
tant  role. 

Our  goal  is  a  unified  model  for  both  earthquake  and 
nuclear  test  data.  The  log  PILg  ratios  are  systematically 
lower  for  earthquakes  than  for  explosions,  which  is  the 
basis  of  the  log  P/Lg  discriminant.  In  order  to  combine 


the  data  sets,  we  subtract  an  average  value  from  each  of 
the  three  data  sets  to  emphasize  the  path  effect.  The 
explosion  data  involve  16  paths,  and  the  earthquake  data 
involve  34  paths.  We  regress  the  three  data  sets  together 
with  n  =  6  explanatory  variables  and  m  =  50  observa¬ 
tions,  16  of  which  are  path  averaged. 

Table  la  lists  the  correlation  matrix  of  log  PILg  ratios 
(Y)  and  six  parameters  that  we  choose  to  characterize 
the  crust  and  Table  lb  lists  the  result  of  multivariate 
regression.  The  definitions  of  the  crustal  parameters  are 
x(\ ),  minimum  altitude;  x(2),  mean  altitude;  x(3), 
minimum  crustal  thickness;  *(4),  mean  crustal  thickness; 
x(5),  maximum  sediment  thickness;  and  a  (6),  yA,  the 
attenuation  term. 

Under  the  assumption  that  the  log  PiLg  ratios  are  a 
linear  function  of  the  six  explanatory  variables  a(1)  to 
x  (6),  we  postulate  the  model 

f  =  B0  +  j^Bj  5*0  )  +  t 
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MINIMUM  Lg  CODA  Q  VALUE 
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Figure  7.  The  log  P/Lg  ratio  variations  with  (a)  maximum  and  (b)  mean  sediment  thickness,  and  (c) 
minimum  Q  value  and  (d)  the  attenuation  term  yA  for  earthquake  data. 


where  j  =  and  t  are  all  vectors  with  50 

components.  Bj,  j  =  1,...,6  are  the  regression 

coefficients  and  Bq  is  the  intercept  term.  The  random 
errors  eiy  i  =  1,...,50  are  assumed  to  be  mutually 
independent  and  to  follow  a  normal  distribution  with 
zero  mean.  In  the  following  steps,  we  examined  et  and 
found  they  are  in  fact  close  to  random  variables  with 
zero  mean. 

The  correlation  between  minimum  altitude  and  mean 
altitude,  *(1)  and  x(2),  is  0.939,  and  that  between 
minimum  and  mean  crustal  thickness,  x(3)  and  *(4),  is 
0.887  (Table  la).  This  means  that  the  local  path  proper¬ 
ties  and  mean  path  properties  are  coupled  in  our  data 
sets.  High  correlations  among  these  variables  imply  that 
multi  col  linearity  exists,  which  may  cause  false  parameter 
estimation  and  instability  of  the  model.  The  last  row  in 
Table  la  gives  the  correlations  of  log  PILg  ratios  (Y) 
with  each  factor,  essentially  combining  the  data  from 
Figures  6  to  11.  The  last  two  factors,  x (5),  maximum 
sediment  thickness  and  *(6),  yA,  have  the  highest  corre¬ 


lations  with  Y .  This  anticipates  the  final  model  we  obtain 
later. 

A  measure  of  the  overall  collinearity  of  independent 
variables  in  multiple  regression  analysis  can  be  obtained 
by  computing  the  condition  number  of  the  correlation 
matrix.  The  condition  number  is  defined  as  the  square 
root  of  the  ratio  of  the  maximum  eigenvalue  to  the 
minimum  eigenvalue  of  the  correlation  matrix.  A  large 
condition  number  provides  evidence  for  strong  collinear¬ 
ity.  The  onset  of  instability  in  a  regression  model  due  to 
collinearity  has  been  empirically  determined  to  occur 
when  condition  numbers  exceed  15  [Chatterjee  and 
Price ,  1991],  In  our  case,  the  condition  number  of  the 
con-elation  matrix  is  found  to  be  10.782.  This  large  con¬ 
dition  number  supports  our  contention  of  collinearity; 
however,  it  is  not  so  severe  as  to  adversely  affect  our 
analysis.  Another  index  to  measure  collinearity  is  the 
sum  of  the  reciprocals  of  the  eigenvalues  (SRE).  As  a 
rule  of  thumb,  if  SRE  is  greater  than  5  times  the  number 
of  explanatory  variables  in  the  problem,  the  variables  are 
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Figure  8.  The  log  P I Lg  ratio  variations  with  (a)  minimum  and  (b)  mean  altitude,  (c)  minimum  and 
(d)  mean  crustal  thickness  for  Novaya  Zemlya  explosion  data. 


collinear  [Chatterjee  and  Price ,  1991].  In  our  case  the 
SRE  value  of  50.065  is  larger  than  30  (5x6),  again 
confirming  the  existence  of  variable  collinearity.  The 
presence  of  collinearity  requires  us  to  assess  the  contri¬ 
bution  of  each  independent  variable  to  the  regression  to 
establish  its  importance.  We  proceed  with  a  backward 
elimination  method  that  begins  with  all  independent  vari¬ 
ables  and  successively  drops  one  at  a  time,  based  on 
their  redundancy. 

The  regression  results  for  all  six  parameters  are  shown 
in  Table  lb,  where  Bj  are  the  regression  coefficients  for 
each  variable  and  S.E.  are  the  standard  errors  associated 
with  each  coefficient.  The  coefficient  of  determination 
R 2  is  used  to  assess  the  fit.  Its  definition  is 

m 

1(37  -J,)2 

R2  =  1  -  — - 

m 

£0 v  -Ji)2 

/  =  ) 

where  y(  are  predicted  values  of  the  dependent  variable, 
yj  is  the  average  value  and  m  is  the  number  of 


observations.  This  quantity  can  be  interpreted  as  the 
proportion  of  the  variability  of  Y  explained  by  the 
regression  on  variables  jc(1  (6).  In  Table  lb, 
R2  =  0.387,  which  says  that  38.7%  of  the  variance  of  Y 
is  accounted  for  by  the  variables  *(1)  to  x(6).  We  do 
not  want  the  value  of  R2  to  significantly  decrease  as  we 
successively  eliminate  variables.  Another  quantity  R2 
(adjusted./?  squared)  is  also  used  [Chatterjee  and  Price , 
1991]  to  assess  goodness  of  fit: 

m 

!>/  1) 

R2  =  \  -  —m - 

£C v,-  -  y,  )2  /  (m— i) 

1=1 

where  n  is  the  number  of  independent  variables.  This 
parameter  is  the  coefficient  of  determination  adjusted  to 
account  for  the  number  of  variables  in  the  regression 
relation.  We  use  both  R2  and  R2  to  provide  stopping 
criteria  in  our  backward  elimination  of  the  variables. 
Here  R2  =  0.031,  which  is  very  low,  and  we  want  to 
maximize  its  value  as  we  eliminate  parameters. 
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Figure  9.  The  log  PILg  ratio  variations  with  (a)  maximum  and  (b)  mean  sediment  thickness,  (c) 
minimum  Q  value  and  (d)  attenuation  term  yA  for  Novaya  Zemlya  explosion  data. 


To  assess  the  significance  of  the  model  that  accounts 
for  38.7%  of  the  variation,  we  use  the  overall  F  value, 
which  has  the  F  distribution  with  n  and  (m-n-1) 
degrees  of  freedom,  as  the  test  statistic  [Chatterjee  and 
Price ,  1991].  Here  n  is  the  number  of  explanatory  vari¬ 
ables  and  m  is  the  number  of  observations.  The  null 
hypothesis  is  all  the  coefficients  Bj,j  =  1,  6  are  zero. 
This  implies  that  there  is  no  linear  relationship  between 
log  P/Lg  ratios  and  any  of  the  path  variables.  The  value 
of  overall  F  computed  from  the  six  variable  regression 
model  is  4.5.  The  1%  right-tail  point  of  the  F  distribu¬ 
tion  for  6  and  43  degrees  of  freedom  is  3.25.  Since  the 
observed  F  value  is  larger  than  this,  the  null  hypothesis 
is  rejected  at  the  99%  significance  level:  all  the 
coefficients  cannot  be  taken  as  zero.  There  is  a  linear 
relationship  between  log  P/Lg  ratio  and  at  least  one  of 
the  path  variables. 

To  determine  which  of  the  independent  variables  is 
the  most  significant  in  the  regression,  partial  F  tests  are 
performed  with  1  and  m-n- 1  degrees  of  freedom.  In  the 
case  of  the  partial  F  test,  the  null  hypothesis  is  that  one 


coefficient,  Bj,  equals  zero  and  that  a  linear  relationship 
exists  between  the  dependent  variable  and  the  remaining 
n- 1  independent  variables.  The  last  column  of  Table  lb 
lists  values  for  the  partial-F  statistic  F(Bj  ~  0).  Only 
the  value  for  x(6)  requires  the  null  hypothesis  to  be 
rejected  at  the  95%  confidence  level,  because  F(B§  =  0) 
is  5.023,  which  is  larger  than  4.07,  the  5%  right-tail 
point  of  the  F  distribution  for  1  and  43  degrees  of  free¬ 
dom.  Therefore  any  of  the  coefficients,  except  x(6)  can 
be  zero.  The  F  statistic  of  x(2)  is  the  lowest  (0.002),  so 
we  decide  to  eliminate  x(2),  mean  altitude,  to  reduce  the 
collinearity  in  the  first  step.  The  variable  to  be  elim¬ 
inated  is  marked  with  ^MPF  (minimum  partial  F)  in 
Table  lb,  and  subsequent  tables. 

After  deleting  x(2),  we  perform  the  regression  again 
with  the  five  remaining  variables.  The  result  is  in  Table 
2.  Now  the  collinearity  is  substantially  reduced.  The 
condition  number  (7.40)  is  32%  less  than  before,  and  the 
sum  of  the  reciprocals  of  eigenvalues  (SRE  =  25.078)  is 
reduced  50.4%,  however,  the  value  is  still  greater  than 
25,  five  times  the  number  of  variables.  Having  removed 
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Figure  10.  The  log  PILg  ratio  variations  with  (a)  minimum  and  (b)  mean  altitude,  (c)  minimum  and 
(d)  mean  crustal  thickness  for  Semipalatinsk  explosion  data. 


variable  x(2 ),  R 2  does  not  decrease,  which  means  the  fit 
remains  almost  the  same.  Meanwhile,  R 2  increases  to 
0.053,  which  represents  improvement  of  the  regression. 
Now  x(l)  has  the  lowest  F  value,  0.024.  Omitting  it 
will  have  the  least  influence  on  the  fit.  We  proceed  to 
do  so  to  seek  a  better  model.  The  resulting  regression  is 
listed  in  Table  3.  The  corresponding  condition  number  is 
6.349  and  SRE  =  18.60. 

Given  the  result  shown  in  Table  3,  we  proceed  to 
remove  jc(3),  and  continue  the  backward  elimination  pro¬ 
cedure.  The  next  result  is  in  Table  4,  The  condition 
number  (2.177)  and  SRE  (4.496)  are  reasonably  small 
this  time.  However,  we  need  to  check  whether  the  esti¬ 
mates  of  all  the  coefficients  are  statistically  robust.  Two 
of  the  variables  still  have  partial  F  statistics  less  than 
4.07,  the  5%  right-tail  point  of  F(l,46),  indicating  that 
another  variable  can  be  eliminated.  Our  two  criteria  for 
stopping  backward  elimination,  a  significantly  decreasing 
R2  and  falling  R2,  have  not  been  obtained,  allowing  us 
to  continue  the  process. 

Table  5  is  the  result  for  only  two  variables,  a  (5),  the 
maximum  sediment  thickness,  and  x(6),  the  yA  term. 


The  condition  number  (1.820)  and  SRE  (2.807)  are  both 
small  now.  The  most  important  result  is  that  both  the 
partial  F  values  of  x(5)  and  x(6)  are  beyond  the  5% 
right-tail  point  of  F(l,47)  =  4.05.  The  estimates  for  these 
two  coefficients  are  now  significant.  This  is  our  pre¬ 
ferred  model.  Its  R2  is  only  4.7%  smaller  than  R2  of  the 
first  model,  and  R2  reaches  its  peak. 

We  did  consider  models  for  single  variables,  x(5)  or 
x (6).  However,  both  the  R2  (the  triangles  in  Figure  12), 
and  R2  (the  solid  circles)  decrease  abruptly.  A  single 
variable  case,  using  either  maximum  sediment  thickness 
or  yA,  does  not  fit  as  well  as  the  two  variable  model. 
This  result  reaffirms  the  importance  of  both  crustal 
attenuation  and  sedimentary  basin  controls  on  regional 
phases  such  as  Lg.  However,  we  emphasize  that  the 
other  crustal  parameters  undoubtedly  play  some  role,  and 
our  statistical  test  merely  identifies  the  most  significant 
factors. 

Application  for  Seismic  Discrimination 

Due  to  the  relatively  larger  S  wave  energy  excited  by 
the  source,  the  Lg  amplitudes  measured  from  earthquake* 
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recordings  are  usually  larger  than  those  of  explosion 
recordings  with  comparable  P  wave  amplitudes.  The 
P  ILg  ratio  is  thus  a  promising  discriminant  for  use  at 
regional  distances  [e.g.,  Bennett  and  Murphy ,  1986;  Tay¬ 
lor  et  al. ,  1989].  However,  variations  in  crustal  struc¬ 
ture  and  attenuation  cause  large  scatter  of  P  !Lg  ratios, 
and  the  discriminant  does  not  work  uniformly  [e.g., 
Lynnes  and  Baumstark ,  1991].  Figure  13a  shows  the 


Table  la.  Correlation  Coefficients 


Variable 

x(l) 

x(2) 

AT  (3) 

x(4) 

*(5) 

jc(6) 

y 

*(n 

1.000 

*(2) 

0.939 

1. 000 

*(3) 

0.805 

0.784 

1.000 

*(4) 

0.864 

0.908 

0.887 

1.000 

*(5) 

-0.529 

-0.541 

-0.549 

-0.517 

1.000 

x(6) 

-0.396 

-0.292 

-0.590 

-0.350 

0.536 

1.000 

Y 

-0.390 

-0.362 

-0.473 

-0.401 

0.499 

0.558 

1.000 

Condition  number  is  10.782  and  the  sum  of  the  reciprocals 
of  eigenvalues  is  50.065. 


log  P!Lg  ratios  for  earthquakes  (solid  circles),  Semipala- 
tinsk  explosions  (crosses),  and  Novaya  Zemlya  explo¬ 
sions  (pluses)  as  a  function  of  m^.  The  two  groups  of 
points  overlap  significantly.  There  may  be  many  reasons 
for  the  overlap,  including  the  physical  separation  of 
paths  recording  explosion  and  earthquake  signals,  the 
relatively  low  frequencies  at  which  the  measurements  are 
made  (P lLg  discriminants  work  best  above  5  Hz,  outside 


Table  lb.  Six-Parameter  Regression 


Variable 

Bj 

S.E. 

F(Bj=0) 

*a> 

0.016 

0.245 

0.005 

x(2) 

0.008 

0.201 

0.002<— MPF 

*(3) 

0.006 

0.026 

0.058 

*(4) 

-0.020 

0.032 

0.393 

x{5) 

0.014 

0.011 

1.668 

*(6) 

0.081 

0.036 

5.023 

Intercept 

-0.138 

R2  -  0.387  R2  =  0.031,  number  of  observations  is  50,  and 
F  (6, 43)  =  4.5. 


ZHANG  ET  AL.:  WAVEGUIDE  EFFECTS  ON  REGIONAL  DISCRIMINANTS 


21,941 


Table  2.  Five-Parameter  Regression 


Variable 

Bj 

S.E. 

F(B;= 0) 

*(  i  > 

0.024 

0.154 

0.024<— MPF 

x(3) 

0.006 

0.025 

0.058 

*(4) 

-0.019 

0.027 

0.514 

*(5) 

0.014 

0.010 

1.793 

x(  6  ) 

0.082 

0.035 

5.469 

Intercept 

0.120 

R2=  0.387,  R2  -  0.053,  number  of  observations  is  50,  and 
F(5,44)  =  5.5. 


the  passband  of  this  data),  as  well  as  possible  failure  of 
the  discriminant  due  to  spall  from  the  explosions.  We 
explore  whether  our  empirical  path  corrections  can 
improve  the  performance  of  the  discriminant  at  all,  given 
this  most  demanding  set  of  conditions.  We  introduce  a 
parameter,  OV,  to  describe  the  overlap  of  the  two  popu¬ 
lations  of  events.  We  speculate  that  the  explosion  data 
cover  a  range  |ii±3aj,  and  earthquake  data  occupy 
Jl2±3a2,  where  |i  is  the  respective  average  ratio  and  a  is 
the  standard  deviation  for  each  population.  We  assume 
JI2  <  \L\*  The  overlap  parameter  is  defined  as 

OV  =  -  +  3g2  - (^' ' 3g,) 

S 

where 

I  (n i  -  l)Gj2  +  (n2  -  1)g22 
S  “  V  n  x . +  *  2 . -  2 

is  the  pooled  standard  deviation  of  the  two  samples 
under  consideration.  Thus  the  overlap  is  expressed  in 
units  of  the  standard  deviation  [Flury  and  Riedwyl , 
1988].  The  n{  and  n2  denote  the  sample  size.  If  the  two 
samples  have  no  overlap  in  range,  the  OV  value  is  nega¬ 
tive.  OV  is  2.64  in  Figure  13a,  which  reflects  significant 
overlap  and  poor  performance  of  the  discriminant. 

Averaging  ratios  for  recordings  from  the  same  source 
can  partially  suppress  path  effects.  When  this  is  done, 
OV  =  0.25,  as  shown  in  Figure  13b.  This  includes 
event-averaged  values  for  all  17  earthquakes  and  97 
explosions.  The  reduction  in  overlap  achieved  by  event 
averaging  indicates  the  extent  of  path  variability, 
although  discrimination  is  still  not  achieved  in  this  case. 

We  use  our  final  model  with  two  explanatory  vari¬ 
ables,  maximum  sediment  thickness  and  yA,  given  in 
Table  5  to  correct  the  log  P/Lg  ratios.  The  reference 


Table  3.  Four-Parameter  Regression 


Variable 

Bj 

S.E. 

F{Bj= 0) 

x(3) 

0.006 

0.025 

0.064<— MPF 

x(  4  ) 

-0.017 

0.023 

0.560 

*(5) 

0.014 

0.010 

1.808 

•*(  6  ) 

0.081 

0.034 

5.566 

Intercept 

0.025 

R2  ~  0.386,  R2  -  0.074,  number  of  observations  is  50,  and 
F  (4, 45)  =  7.1. 


Table  4.  Three-Parameter  Regression 


Variable 

Bj 

S.E. 

o 

li 

sT 

if 

x(4) 

-0.012 

0.010 

1.31  lf-MPF 

x(5  ) 

0.014 

0.010 

1.920 

*(6) 

0.076 

0.027 

8.119 

Intercept 

0.000 

R2  =  0.386,  R2  ~  0.093,  number  of  observations  is  50,  and 
F(3,46)  =  9.6. 


sediment  thickness  is  chosen  as  10  km,  and  the  reference 
yA  is  4.  The  choice  of  reference  values  is  arbitrary,  and 
does  not  affect  the  result.  If  h  denotes  sediment  thick¬ 
ness,  the  correction  is  done  by: 

Ratiocorrected  =  Ratioraw  -  (0.0195 h  +  0.0795(yA)). 

After  correction  for  the  path  effect,  the  variance  of  the 
ratio  values  for  earthquake  data  is  reduced  from  0.3078 
to  0.1846,  about  40%.  The  variance  for  explosion  obser¬ 
vations  is  reduced  from  0.1648  to  0.1202,  more  than 
27%.  The  OV  value  decreases  from  2.64  to  1.47.  This 
change  is  still  not  enough  to  separate  the  two  groups,  but 
it  is  a  significant  improvement.  The  event  average  values 
in  Figure  lid  also  show  reduced  overlap.  We  do  not 
explore  more  quantitative  discriminant  measures  with 
this  data  because  it  is  clear  that  even  with  the  path 
corrections  we  have  not  eliminated  overlap.  It  will  be  of 
great  interest  to  establish  whether  comparable  path 
correction  procedures  improve  discriminant  performance 
for  higher  frequency  P/Lg  data  and/or  for  explosion  and 
earthquake  data  sampling  similar  paths. 

Discussion  and  Conclusions 

The  empirical  relationships  established  in  this  paper 
and  by  Zhang  and  Lay  [1994a,b]  indicate  that  regional 
wave  propagation  is  influenced  by  attributes  of  gross  cru¬ 
stal  structure  that  can  be  independently  determined.  This 
has  important  implications  for  both  quantifying  the  wave 
propagation  effects,  as  well  as  for  developing  predictive 
ability  for  regional  phase  energy  partitioning  on  new 
paths  through  the  crust.  The  latter  ability  is  of  particular 
importance  for  discriminating  explosion  and  earthquake 
sources  in  regions  with  no  prior  known  explosion 
activity.  At  this  time,  numerical  methods  such  as  finite 
difference  are  inadequate  to  perform  extensive  modeling 
of  the  heterogeneous  waveguide  effects  indicated  by  the 
crustal  models  used  here,  but  the  trends  found  above  will 


Table  5.  Two-Parameter  Regression 


Variable 

Bj 

S.E. 

o 

it 

2? 

x(5) 

0.019 

0.009 

4.172 

x(6) 

0.079 

0.027 

8.823 

Intercept 

-0.568 

R2  -  0.368,  R2  =  0.099,  number  of  observations  is  50,  and 
F(2,47)  =  13.7. 
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Figure  12.  The  behavior  of  the  coefficient  of  determination  R2  and  adjusted  R  -square  R2  as  func- 
tions  of  backward  elimination  step.  The  triangles  indicate  R2,  which  drops  only  a  tiny  amount  in 
t  e  first  four  steps,  but  then  drops  strongly  in  the  fifth  step,  when  we  use  only  one  explanatory  vari- 
able,  sediment  depth  or  attenuation.  The  solid  circles  are  R2,  which  rises  until  the  fourth  variable 
elimination  and  then  drops.  These  two  statistics  provide  the  stopping  criteria  of  backward  elimina- 
tion  for  multivariate  regression. 


guide  future  modeling  efforts,  as  well  as  provide  empiri¬ 
cal  corrections. 

The  multivariate  analysis  supports  the  importance  of 
both  waveguide  heterogeneity  and  crustal  attenuation  in 
shaping  the  regional  signal.  Use  of  higher-resolution  cru¬ 
stal  models  and  more  accurate  attenuation  models  may 
reveal  the  important  physics  better  than  this  initial  study. 
In  addition,  most  of  our  propagation  paths  are  long,  and 
there  is  little  crossing  coverage  with  which  to  isolate 
effects. 

We  use  the  Lg  coda  Q  model  of  Pan  et  al  [1992]  as 
a  parameterization  of  attenuation  effects  rather  than  as 
an  explicit  Q  model.  The  Lg  coda  Q  model  undoubtedly 
reflects  some  of  the  Lg  Q  effects  on  any  given  path,  but 
it  is  likely  that  it  will  not  provide  a  full  characterization 
of  the  apparent  attenuation  for  individual  paths.  In 
deriving  the  Lg  coda  Q  model,  the  Lg  coda  waves  are 
assumed  to  be  Lg  energy  scattered  within  an  ellipse,  fol¬ 
lowing  Xie  and  Nuttli  [1988].  Thus  the  Q  value  should 
reflect  only  the  attenuation  properties  of  the  crust.  How¬ 
ever,  there  may  be  some  mantle  contribution.  Figure  14 
shows  empirical  relationships  between  the  attenuation 


term  yA  and  four  phases  from  the  earthquake  data:  RMS 
P  amplitudes  measured  in  the  time  window  from  0  to  20 
s,  P  coda  from  20  to  50  s,  Sn ,  and  Lg .  The  log  Lg  is 
strongly  correlated  with  yA  (Figure  14d),  but  the  absolute 
value  of  the  slope  (-0.27)  is  much  smaller  than  the 
theoretical  value  for  1  Hz  of  \og]0e  =  -0.43.  This  indi¬ 
cates  that  the  exponential  attenuation  relationship  for  a 
single  frequency  (1  Hz)  cannot  be  directly  used  for  the 
band-passed  data  (0.6-3  Hz  in  our  study).  The  correla¬ 
tion  of  yA  with  log  Sn  is  the  second  strongest.  Sn  waves 
mainly  propagate  in  the  upper  mantle  but  do  involve  cru¬ 
stal  paths.  This  correlation  with  the  attenuation  term 
implies  that  the  Lg  coda  used  to  determine  the  Q  model 
may  sample  the  upper  mantle  as  well.  However,  Sn  to 
Lg  scattering  may  be  responsible  for  the  Sn  coda  having 
sensitivity  to  crustal  parameters.  [Baumgardt ,  1990]. 
Although  the  direct  P  waves  (0-20  s,  Figure  14a)  corre¬ 
late  weakly  with  yA,  the  P  wave  coda  (20-50  s,  Figure 
14b)  correlation  is  surprisingly  strong.  Perhaps  scattered 
Lg  contributes  to  the  P  wave  coda,  or  again  there  may 
be  some  mantle  contribution  to  the  Q  values.  The  corre¬ 
lation  of  P  (0-50  s),  used  above,  is  almost  as  strong  as 


ZHANG  ET  AL.:  WAVEGUIDE  EFFECTS  ON  REGIONAL  DISCRIMINANTS 


21,943 


Figure  13.  The  log  PILg  ratios  for  earthquakes  (solid  circles)  and  explosions  (crosses  for  Semipala- 
tinsk  data  and  pluses  for  Novaya  Zemlya  data)  as  a  function  of  m^.  (a)  individual  raw  data,  (b) 
event-averaged  raw  data,  (c)  individual  data  after  applying  the  multivariate  model  to  correct  for  path 
effects,  and  (d)  event-averaged  corrected  data.  The  variances  and  overlaps  of  both  groups  are 
reduced  after  correction. 


P  (20-50  s).  It  appears  that  the  Q  model  generally  has 
the  right  geographic  pattern,  and  certainly  attenuation  is 
an  important  factor  influencing  the  regional  wave  propa¬ 
gation  for  the  various  phases.  However,  it  seems  that  Lg 
coda  Q  not  only  represents  the  attenuation  property  of 
crust  but  also  that  of  upper  mantle,  given  the  strong  rela¬ 
tions  with  Sn  and  P  coda.  This  may  also  indicate  cou¬ 
pling  between  Q  structure  in  the  crust  and  upper  mantle, 
a  natural  consequence  of  tectonic  processes.  Unlike  for 
Lg ,  we  do  not  have  any  theoretical  framework  for 
directly  correcting  the  P  amplitudes  for  attenuation,  thus 
we  have  treated  the  effect  of  yA  on  the  PiLg  values 
empirically. 

Zhang  and  Lay  [1994b]  attribute  correlations  found 
with  maximum  water  depth  to  either  locally  thick  sedi¬ 
ment  layers  or  thin  crust  and  could  not  find  a  robust  pro¬ 
cedure  to  separate  the  effects  of  the  two  factors.  One  of 
the  objectives  of  this  study  was  to  better  elucidate  propa¬ 
gation  effects  on  regional  signals.  We  obtain  a /'best  fit" 
model  through  multivariate  regression;  however,  we 


emphasize  that  the  solution  found  is  not  unique.  Besides 
waveguide  attenuation,  the  structural  parameters  may  all 
have  some  importance.  The  maximum  sediment  thick¬ 
ness  factor  (x(5))  is  correlated  with  the  other  structural 
factors  with  correlation  coefficients  >  0.5  (Table  la).  If 
we  use  yA  plus  any  of  the  other  structural  factors,  instead 
of  maximum  sediment  thickness,  the  resulting  fits  are  not 
far  from  the  best  one.  Table  6  shows  these  results,  where 
the  columns  indicate  models  with  yA  plus  one  of  the  fol¬ 
lowing:  minimum  altitude,  mean  altitude,  minimum  cru¬ 
stal  thickness,  mean  crustal  thickness,  and  maximum 
sediment  thickness.  All  the  coefficients  of  determination 
R 2  are  greater  than  0.344,  only  slightly  smaller  than 
0.368  for  the  best  model  with  maximum  sediment  thick¬ 
ness  and  yA.  All  the  overall  F  values  are  greater  than 
the  1%  right-tail  point  of  F{ 2,  46)  =  5.10.  However, 
only  the  maximum  sediment  thickness  coefficient  passes 
the  partial  F  test;  in  the  row  marked  "partial  F'\  4.172 
>  5%  right-tail  point  of  F(l,  46)  =  4.05.  Our  choice  of 
the  "best  model"  with  yA  and  maximum  sediment  thick- 
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Figure  14.  Correlations  of  the  path  attenuation  term  yA,  computed  from  the  Lg  coda  Q  model 
[Pan  et  al ,  1992]  with  (a)  P  wave  RMS  amplitude  measured  in  the  0-20  s  window,  (b)  P  wave 
coda  RMS  amplitude  measured  in  a  20-50  s  window,  (c)  RMS  Sn ,  and  (d)  RMS  Lg .  These  correla¬ 
tions  imply  that  the  Lg  coda  Q  reflects  not  only  crustal  properties  but  also  the  upper  mantle  to  a 
lesser  extent. 


ness  is  clearly  preferred  in  a  statistical  sense  but  the 
strong  collinearity  in  the  problem  should  not  be  over¬ 
looked. 

Although  the  crustal  thickness  is  not  included  in  our 
final  model,  its  importance  cannot  be  neglected.  Its  par¬ 
tial  F  value  is  only  30%  less  than  that  of  sediment 
thickness  (Table  4).  Our  study  area  includes  two  regions 
with  some  of  the  thickest  sediment  layers  in  the  world, 
the  effect  of  which  is  pronounced.  In  other  regions  lack¬ 
ing  such  thick  sediment  regions,  crustal  thickness  meas- 


Table  6.  Comparison  of  Alternate  Models 


yA  + 

minimum 

altitude 

mean 

altitude 

minimum 

crustal 

thickness 

mean 

crustal 

thickness 

maximum 

sediment 

thickness 

R 2  0.346 

0.355 

0.344 

0.360 

0.368 

Over  all  F  12.4 

13.0 

12.3 

13.2 

13.7 

Partial  F  2.428 

3.167 

2.273 

3.521 

4.172 

ures  may  be  more  influential.  Measures  of  crustal  thick¬ 
ness  variation  should  be  further  explored  for  their 
influence  on  the  regional  signals. 

While  we  gain  some  insight  into  regional  wave  energy 
partitioning  by  these  empirical  approaches,  synthetic  cal¬ 
culations  are  ultimately  needed  to  fully  understand  the 
coupled  relationships.  The  coupling  of  structural  parame¬ 
ters  and  Q  is  also  important  to  consider.  Mitchell  and 
Hwang  [1987]  conclude  that  Lg  Q  and  sediment  thick¬ 
ness  are  coupled  in  the  central  United  States,  but  are  not 
related  in  the  western  United  States.  It  is  possible  that 
these  factors  are  related  in  portions  of  Eurasia  and  not  in 
other  regions,  giving  rise  to  significance  as  two  separate 
variables  in  our  multivariate  regressions.  Variation  in 
sediment  thickness  may  also  play  an  important  role  in  Lg 
blockage. 

Our  primary  conclusion  is  that  the  combined  effect  of 
the  two  crustal  variables:  yA,  which  reflects  the  gross 
attenuation  properties  of  each  path,  and  maximum  sedi¬ 
ment  thickness,  which  is  the  most  prominent  structure 
factor  causing  blockage  of  Lgy  dominate  the  path  effect 
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on  the  regional  log  PILg  ratio.  This  offers  the  potential 
to  reduce  variance  in  measurements  of  these  ratios  for 
seismic  discrimination  applications. 
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